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ABSTRACT 

This  report  forms  the  user’s  guide  for  Version  1.0  of  LSSOL,  a  set  of  Fortran  77  subroutines  for 
linearly  constrained  linear  least-squares  and  convex  quadratic  programming.  The  method  of  LSSOL 
is  of  the  two-phase,  active-set  type,  and  is  related  to  the  method  used  in  the  package  SOL/QPSOL « 
(diU  et  al.,  1984b^-Two  main  features  of  LSSOL  are  its  exploitation  of  convexity  and  treatment 
of  singularity. 

LSSOL  may  also  be  used  for  linear  programming,  and  to  find  a  feasible  point  with  respect  to  a 
set  of  linear  inequality  constraints.  LSSOL  treats  all  matrices  as  dense,  and  hence  is  not  intended 


for  large  sparse  problems. 


t  LSSOL  is  available  from  the  Stanford  Office  of  Technology  Licensing,  350  Cambridge  Avenue, 
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1.  PURPOSE 

LSSOL  is  ii  collection  of  Fortran  77  subroutines  designed  to  solve  a  class  of  quadratic  prograniiniiig 
problems  that  are  assumed  to  be  stated  in  the  following  general  form: 


LCLS 

minimize 

F{x) 

subject  to 

where  C  is  mj. 

X  n  (m^  may  be  zero)  and  F{x)  is 

one  of  the  following  objective  functions: 

FP: 

None 

(find  a  feasible  point  for  the  constraints) 

LP: 

c^x 

QPl: 

\x'^Ax 

A  symmetric  and  positive  semi-definite, 

QP2: 

c^x  -f-  \x^Ax 

A  symmetric  and  positive  semi-definite, 

QP3: 

^x^A^Ax 

A  m  X  n  upper-trapezoidal. 

qP4: 

c^x  -t-  ^x^A'^Ax 

A  m  X  n  upper-trapezoidal. 

LSt: 

illb-Axll^ 

A  m  X  n, 

LS2; 

c^x  -h  j||b  -  Ax\\^ 

A  m  X  n, 

LS3: 

A  m  X  n  upper-trapezoidal, 

LS4: 

|ll^-  ^*11* 

A  m  X  n  upper- trapezoidal, 

with  c  an  n-vcctor  and  b  an  m-vector.  The  specific  objective  function  to  be  minimized  is  selected 
using  the  optional  parameter  Problem  Type  (see  Section  4.2).  In  all  that  follows,  problems  of 
type  “LP”,  “QP’’  and  “LS”  will  be  referred  to  as  linear  programming,  quadratic  programming  and 
constrained  least-squares  problems  respectively. 

The  constraints  involving  C  will  be  called  the  general  constraints.  Note  that  upper  and  lower 
bounds  are  specified  for  all  the  variables  and  for  all  the  general  constraints.  An  equality  constraint 
is  specified  by  setting  =  Wi.  If  certain  bounds  are  not  present,  the  associated  elements  of  f  or  u 
can  be  set  to  special  values  that  will  be  treated  as  — oo  or  -l-oo.  (See  the  description  of  the  optional 
parameter  Infinite  Bound  in  Section  4.2.) 

The  constant  second-derivative  matrix  of  F(x)  is  defined  as  ff,  the  Hessian  matrix.  In  the 
LP  case,  H  —  0.  In  QP  cases  1  and  2,  H  =  A',  and  in  QP  cases  3  and  4,  //  =  A^A.  In  all  LS 
cases,  II  =  A^A.  Problems  of  type  QP3  or  QP4  with  A  not  in  trapezoidal  form  should  be  solved 
as  type  LSI  or  LS2  with  b  =  0.  When  considering  problems  of  type  LS,  we  shall  refer  to  A  as  the 
least-squares  matrix  and  to  b  as  the  vector  of  observations. 

The  user  must  supply  an  initial  estimate  of  the  solution.  If  the  Hessian  matrix  is  non-singular, 
LSSOL  will  obtain  the  unique  (global)  minimum.  If  H  is  singular,  the  solution  may  still  be  a  global 
minitrium  if  all  active  constraints  have  nonzero  Lagrange  nniltipliers.  Otherwise,  the  solution 
obtained  will  either  be  a  weak  minimum  (i.e.,  with  a  nni<pu‘  optimal  objective  value,  but  an 
infinite  set  of  optimal  i),  or  else  the  objective  function  is  unbounded  below  in  the  feasible  region. 
The  last  case  can  occur  only  when  F(x)  contains  an  cxpUcit  linear  term  (as  in  problems  of  type 
LP,  QP2,  QP4,  LS2  and  LS4). 

The  LSSOL  package  contains  approximately  6000  lines  of  ANSI  Fortran  77,  of  whicli  about 
50%  arc  comments. 
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2.  DESCRIPTION  OF  THE  ALGORITHM 

Here  we  briefly  smiiui.-iri/.e  the  inaiii  fciilures  of  the  method  (if  LSSOL.  Wliere  possible,  explicit 
reference  is  made  to  the  names  of  variables  that  are  parameters  of  subroutine  LSSOL  or  appear  in 
the  printed  output. 

The  method  of  LSSOL  is  a  two-pha.se  (primal)  quadratic  programming  method  (see  Gill  ct  hL, 
198  lb)  with  features  to  exploit  the  convexity  of  the  objective  function.  (In  the  full-rank  case,  the 
method  is  related  to  that  of  Stoer,  1971.)  The  two  phases  of  the  method  are:  flndiiig  an  initial 
feasible  point  by  minimizing  the  sum  of  infeiisibilities  (the  feasibility  phase),  and  minimizing  the 
(piadratic  objective  function  within  the  feasible  region  (the  optimality  phase).  The  computations 
in  both  phases  are  performed  by  the  same  subroutines.  The  two-phase  nature  of  the  algorithm  is 
reflected  by  changing  the  function  being  minimized  from  the  sum  of  infeasibilities  to  the  quadratic 
objective  function.  The  feasibility  phase  does  not  perform  the  standard  simplex  method  (i.e.,  it 
does  not  necessarily  find  a  vertex),  except  in  the  LP  case  when  <  n.  Once  any  iterate  is  feasible, 
all  subsequent  iterates  remain  feasible. 

In  general,  an  iterative  process  is  required  to  solve  a  quadratic  program.  (For  simplicity,  we 
shall  always  consider  a  typical  iteration  and  avoid  reference  to  the  index  of  the  iteration.)  Each 
new  iterate  x  is  defined  by 

X  =  x  +  ap,  (1) 

where  the  step  length  a  is  a  non-negative  scalar,  and  p  is  called  the  search  direction. 

At  each  point  x,  a  working  set  of  constraints  is  defined  to  be  a  linearly  independent  subset 
of  the  constraints  that  arc  satisfied  “exactly”  (to  within  the  tolerance  defined  by  the  optional 
parameter  “Feasibility  Tolerance”;  see  Section  4.2).  The  working  set  is  the  current  prediction 
of  the  constraints  that  hold  with  equality  at  a  solution  of  LCLS.  The  search  direction  is  constructed 
so  that  the  constraints  in  the  working  set  remain  unaltered  for  any  value  of  the  step  length.  For 
a  bound  constraint  in  the  working  set,  this  property  is  achieved  by  setting  the  corresponding 
component  of  the  search  direction  to  zero.  Thus,  the  associated  variable  is  fixed,  and  specification 
of  the  working  set  induces  a  partition  of  i  into  fixed  and  free  variables.  During  a  given  iteration, 
the  fixed  variables  are  effectively  removed  from  the  problem;  since  the  relevant  components  of  the 
search  direction  are  zero,  the  columns  of  C  corresponding  to  fixed  variables  may  be  ignored. 

Let  m,v  denote  the  number  of  general  constraints  in  the  working  set  and  let  denote  the 
number  of  variables  fixed  at  one  of  their  bounds  (mw  and  Upx  are  the  ejuantities  “Lin”  and  “End” 
in  the  printed  output  from  LSSOL).  Similarly,  let  tIfr  (npR  —  n  —  n^x)  denote  the  number  of  free 
variables.  At  every  iteration,  the  variables  arc  re-ordered  so  that  the  last  variables  are  fixed, 
with  all  other  relevant  vectors  and  matrices  ordered  accordingly.  The  order  of  the  variables  is 
indicate(l  by  the  list  of  indices  KX,  a  parameter  of  LSSOL. 

Let  Cf  H  denote  the  m„.  x  UpR  submatrix  of  general  constraints  in  the  working  set  corresponding 
to  the  free  variables,  and  let  p^R  denote  the  search  direction  with  respect  to  the  free  variables  only. 
The  general  constraints  in  the  working  .set  will  be  unaltered  by  any  move  along  p  if 

C'frPfh  —  0-  (2) 

In  order  to  compute  pp„,  the  TQ  factorization  of  Grr  is  used: 

GprQpr  =  (0  T),  (3) 

where  T  is  a  nonsingular  m„  x  m„  reverse- triang(dar  matrix  (i.e.,  tij  --  0  if  i  -fy  <  m„  ),  and  the 
non-singular  n,R  x  n^R  matrix  Qm  is  the  product  of  orthogonal  transformations  (see  Gill  et  al., 
1981a).  If  the  columns  of  G,  r  arc  partitioned  so  that 

QrR  =  {Z  Y), 


(4) 
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where  Y  is  n^n  x  riiw,  then  the  (n^  =  —  m,,. )  columns  of  Z  form  a  basis  for  the  null  space 

of  CpR.  Thus,  Ppr  will  satisfy  (2)  only  if 


Pfh  —  Zpz 


(5) 


for  some  vector  p^. 

Let  Q  denote  the  n  x  n  matrix 


where  is  the  identity  matrix  of  order  Tipx-  Let  R  denote  an  n  x  n  upper-triangular  matrix  (the 
Chulesky  factor)  stich  that 

-  R^R,  (7) 

and  let  the  matrix  of  first  rows  and  columns  of  R  be  denoted  by  R^.  (Recall  that  R  in  (7)  will 
in  general  have  been  re-ordered.) 

The  definition  of  in  (5)  depends  on  whether  or  not  the  matrix  R^  is  singular  at  x.  In  the 
non-singular  case,  p^  satisfies  the  equations 

=  -9z,  (8) 

where  denotes  the  vector  Z^g^y^  and  g  denotes  the  objective  gradient,  (The  norms  of  is 
the  printed  quantity  Norm  Gf .)  When  p^  is  defined  by  (8),  i  -f  p  is  the  minimizer  of  the  objective 
function  subject  to  the  constraints  (bounds  and  general)  in  the  working  set  treated  as  equalities. 
In  general,  a  vector  is  available  such  that  R^f^  =  —g^,  which  allows  p^  to  be  computed  from 
a  single  back-substitution  R^Pz  =  fz-  For  example,  when  solving  problem  LSI,  comprises  the 
first  elements  of  the  transformed  residual  vector 

f  =  P(b-  Ax),  (9) 

which  is  recurred  from  one  iteration  to  the  next,  where  P  is  an  orthogonal  matrix. 

In  the  singular  case,  p^  is  defined  such  that 

RzPz  =  0  and  <  0.  (10) 

This  vector  has  the  property  that  the  objective  function  is  linear  along  p  and  may  be  reduced  by 
any  step  of  the  form  x  +  ap,  a  >  0. 

The  vector  Z^Pf^  is  known  as  the  projected  gradient  at  x.  If  the  projected  gradient  is  zero, 
X  is  a  constrained  stationary  point  in  the  siibspace  defined  by  Z.  During  the  feasibility  phase,  the 
projected  gradient  will  usually  be  zero  only  at  a  vertex  (although  it  may  be  zero  at  non-vertices  in 
the  presence  of  constraint  dependencies).  Dtiring  the  optimality  phase,  a  zero  projected  gradient 
implies  that  x  minimizes  the  quadratic  objective  when  the  constraints  in  the  working  set  are  treated 
as  equalities,  .^t  a  constrained  stationary  point,  Lagrange  multipliers  Ac  and  A„  for  the  general 
and  bound  constraints  are  defined  from  the  equations 

^FR'^c  UpR  und  A,,  —  —  Crj^Ac-  (H) 

Given  a  positive  constant  6  of  the  order  of  the  machine  precision,  the  Lagrange  multiplier  Ay 
corresponding  to  an  inecjuality  constraint  in  the  working  set  is  said  to  be  optimal  if  Ay  <  A  when 
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the  associated  constraint  is  at  its  upper  hound,  or  if  Xj  >  —6  when  the  associated  constraint  is 
at  its  lower  hound.  If  a  inultijdier  is  non-optinial.  the  objective  function  (either  the  true  r)bjective 
or  the  sum  of  infeasibilities)  can  be  reduced  by  deleting  the  corresponding  constraint  (with  iiulex 
Jdel;  see  Section  5)  from  the  working  set. 

If  optimal  multii)liers  occur  during  the  feasibility  phase  and  the  sum  of  infeasibilities  is  nonzero, 
there  is  no  feasible  point,  and  LSSOL  will  continue  until  the  minimum  value  of  the  sum  of  infeasi¬ 
bilities  has  been  found.  At  this  point,  the  Lagrange  multiplier  Xj  corresponding  to  an  inequality 
constraint  in  the  working  set  will  be  such  that  —(1  +^)  <  Aj  <  ^  when  the  associated  constraint 
is  at  its  upper  hound,  and  —6  <  Aj  <  1  -f-  6  when  the  associated  constraint  is  at  its  lower  hound. 
Lagrange  multipliers  for  e<}unlity  constraints  will  satisfy  |Ajj  <  1-4-6. 

The  choice  of  step  length  is  based  on  remaining  feasible  with  respect  to  the  satisfied  constraints. 
If  Rz  is  nonsingular  and  x  +  p  is  feasible,  a  will  be  taken  as  unity.  In  this  case,  the  projected 
gradient  at  x  will  be  zero,  and  Lagrange  multipliers  are  computed.  Otherwise,  a  is  set  to  the 
step  to  the  “nearest”  constraint  (with  index  Jadd;  see  Section  5),  which  is  added  to  the  working 
set  at  the  next  iteration. 

If  A  is  not  ini)ut  as  a  triangular  matrix,  it  is  overwritten  by  a  triangular  matrix  R  satisfying 
(7)  obtained  using  the  Cholesky  factorization  in  the  QP  case,  or  the  QR  factorization  in  the  LS  case. 
Crdumn  interchanges  are  used  in  both  cases,  and  an  estimate  is  made  of  tlie  rank  of  the  triangular 
factor.  Thereafter,  the  dependent  rows  of  R  are  eliminatecl  from  the  problem. 

Each  change  in  the  working  set  leads  to  a  simple  change  to  Cpr:  if  the  status  of  a  general 
constraint  changes,  a  row  of  Cpi,  is  altered;  if  a  bound  constraint  enters  or  leaves  the  working  set, 
a  co/iinin  of  Cm  changes.  Explicit  representations  are  recurred  of  the  matrices  T,  (?fh  a-nd 

of  vectors  Q^g,  Q^c  and  /,  which  are  related  by  the  formulae 


and 


/  =  e6-(")e^x 


(6  =  0  for  the  QP  case). 


Q^!J  rr.  Q^c  -  RTf, 


Note  that  the  triangular  factor  R  associated  with  the  Hessian  of  the  original  problem  is  updated 
during  both  the  optimality  and  the  feasibility  phases. 


The  treatment  of  the  singular  case  depends  critically  on  the  following  feature  of  the  matrix 
updating  schemes  used  in  LSSOL:  if  a  given  factor  /?,  is  non-singular,  it  can  become  singular 
during  subsecpient  iterations  only  when  a  constraint  leaves  the  working  set,  in  which  case  only  its 
last  diagonal  element  can  become  zero.  This  property  implies  that  a  vector  satisfying  (10)  may 
be  fo\ind  using  the  single  back-substitution  R^p^  =  Cz,  wliere  is  tlie  matrix  R^  with  a  unit 
last  iliagonal,  and  is  a  vector  of  all  zeros  except  in  the  last  position.  If  H  is  singular,  the 
matrix  R  (and  hence  R^)  rnciy  be  singular  at  the  start  of  the  optimality  phase.  However,  R^  will 
be  non-singular  if  enough  constraints  are  included  in  the  initial  working  set.  (The  mdl  matrix  is 
positive  definite  by  definition,  corresponding  to  the  case  when  Cy,,  contains  constraints.)  The 
idea  is  to  include  as  many  general  constraints  as  necessary  to  ensure  a  non-singular  R^. 

At  the  beginning  of  each  phase,  an  upper-triangular  matrix  Ry  is  determined  that  is  the  largest 
non-singular  leading  submatrix  of  R^.  The  use  of  interchanges  during  the  factorization  of  A  tends 
to  maximize  the  dimension  of  /?].  (The  rank  of  R^  is  estimated  using  the  optional  parameter  Rank 
Tolerance;  see  Section  4.2.)  Let  Zj  denote  the  columns  of  Z  corresponding  to  /?i ,  and  let  Z  be 
partitioned  as  Z  =  (  Zj  Z2  ).  A  working  set  for  which  Z^  defines  the  null  space  can  be  obtained 


by  iiidiidin;'  the  rows  of  Zj  'is  “artificial  constraints”.  Minimization  of  the  objective  function  then 
proceeds  within  the  subspace  defined  by  Zi. 

The  artificially  au^'iiiented  working  set  is  given  by 

"  (  z}  )  ’ 

so  that  p,  „  will  satisfy  Cf^Pf^  =  0  and  Z^Pr^  =  0.  By  definition  of  the  TQ  factorization,  Ck,, 
automatically  satisfies  the  following: 

)  (h  n  =  (  ZT  )  (  ^2  ^  )  =  (  0  T  ), 


where 


-c  :)• 


and  h<'nce  the  TQ  factorization  of  (12)  rccpiires  no  additional  work. 

The  matrix  Z2  need  not  be  kept  fixed,  since  its  role  is  purely  to  define  an  appropriate  null  space; 
the  TQ  factorization  can  therefore  be  updated  in  the  normal  fashion  as  the  iterations  proceed. 
No  work  is  required  to  “delete”  the  artificial  constraints  associated  with  when  Z^g^^  =  0, 
since  this  sinqily  involves  repartitioning  Qyn-  When  deciding  which  constraint  to  delete,  the 
“artificial”  multiplier  vector  associated  with  the  rows  of  Z^  is  ecpial  to  Z^OyR,  and  the  multipliers 
cv>rrespondi/ig  to  the  rows  of  the  “true”  working  set  are  the  multipliers  that  wotild  be  obtained  if 
the  temporary  constraints  were  not  present. 

The  number  of  columns  of  Z  and  Zi ,  the  Euclidean  norm  of  Z^g^R ,  and  the  condition  estimator 
of  /?i  appear  in  the  printed  output  as  Nz,  Nzl,  Norm  Gzl  and  Cond  Rzl  (see  Section  5). 

Although  the  algorithm  of  LSSOL  does  not  perform  simplex  steps  in  general,  there  is  one 
exception:  a  linear  program  with  fewer  general  constraints  than  variables  (i.e.,  ni,,  <  n).  (Use 
of  the  simplex  method  in  this  situation  leads  to  savings  in  storage.)  At  the  starting  point,  the 
“natural”  working  set  (the  set  of  constraints  exactly  or  nearly  satisfied  at  the  starting  point) 
is  augmented  with  a  suitable  number  of  “temporary”  boumls,  eacli  of  which  has  the  effect  of 
temporarily  fixing  a  variable  at  its  current  value.  In  subsecpient  iterations,  a  temporary  bound  is 
treated  as  a  standard  constraint  until  it  is  deleted  from  the  working  set,  in  which  case  it  is  never 
added  again. 

One  of  the  most  important  features  of  LSSOL  is  its  control  of  the  conditioning  of  the  working 
set.  wliose  nearness  to  linear  dependence  is  estimated  by  the  ratio  of  the  largest  to  smallest  diagonals 
of  the  TQ  factor  T  (the  printed  value  Cond  T;  see  Section  5).  In  constructing  the  initial  working  set, 
constraints  are  excluded  that  would  result  in  a  large  value  of  Cond  T.  Thereafter,  LSSOL  allows 
constraints  to  be  violated  by  as  much  as  a  user-specified  Feasibility  Tolerance  (see  Section 
4.2)  in  order  to  provide,  whenever  possible,  a  choice  of  constraints  to  be  added  to  the  working  set 
at  a  given  iteration.  Let  denote  the  maximum  step  at  which  x  -f  n^p  does  not  violate  any 
constraint  by  more  than  its  feasibility  tolerance.  All  constraints  at  distance  n  (a  <  along  p 
from  the  current  point  are  then  viewed  as  acceptable  candidates  for  inclusion  in  the  working  set. 
The  constraint  whose  normal  makes  the  largest  angle  with  the  search  direction  is  added  to  the 
working  set.  In  order  to  ensure  that  the  new  iterate  satisfies  the  constraints  in  the  working  set  as 
accurately  as  possible,  the  step  taken  is  the  exact  distance  to  the  newly  added  cr)nstraint.  As  a 
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consequence,  negative  steps  are  occasionally  permitted,  since  the  current  iterate  may  violate  the 
constraint  to  be  added  by  as  much  as  the  feasibibty  tolerance. 

LSSOL  has  been  designed  to  be  efficient  when  used  to  solve  a  sequence  of  related  problems  -  for 
example,  within  a  sequential  quadratic  programming  method  for  nonlinearly  constrained  optimiza¬ 
tion  (e.g.,  the  NPSOL  package  of  Gill  et  a].,  1986).  In  particular,  the  user  may  specify  an  initial 
working  set  (the  indices  of  the  constraints  believed  to  be  .satisfied  exactly  at  the  solution);  see  the 
discussion  of  the  optional  parameter  Warm  Start  in  Section  4.2. 
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3.  SPECIFICATION  OF  SUBROUTINE  LSSOL 

Tlu’  forin:il  specification  of  LSSOL  is  the  following; 

SUBROUTINE  LSSOL  (  H,  N, 

NCLIN.  NROWC,  NROWA, 

C.  BL,  BU.  CVEC, 

ISTATE,  KX,  X.  A.  B. 

INFORM,  ITER.  OBJ,  CLAMDA, 

IW,  LENIW,  W.  LENW  ) 

M.  N,  NCLIN 

NROWG,  NROWA,  INFORM,  ITER,  LENIW,  LENW 
ISTATE(N+NCLIN).  KX(N),  IW(LENIW) 

OBJ 

CCNROWC,*),  BL(N+NCLIN).  BU(N+NCLIN) , 

CVEC(*),  X(N).  AfNROWA,*), 

B(*),  CLAMDA (N+NCLIN),  W(LENW) 

Note;  Here  and  elsewhere,  the  specification  of  a  parameter  as  REAL  should  be  interpreted  as  working 

prreision,  which  may  be  DOUBLE  in  some  installations. 

3.1.  Formal  parameters 

M  (Input)  The  number  of  rows  in  the  array  A.  If  the  problem  is  specified  as  type  FP  or 

LP  (see  Section  4),  M  is  not  referenced  and  is  assumed  to  be  zero. 

If  the  problem  is  of  type  QP,  M  will  usually  be  N,  the  number  of  variables.  However,  a 
value  of  M  loss  than  N  is  appropriate  for  QP3  or  QP4  if  A  is  an  upper-trapezoidal  matrix 
with  M  rows.  Similarly,  M  may  be  used  to  define  the  dimension  of  a  leading  block  of 
non-zeros  in  the  Hessian  matrices  of  QPl  or  QP2,  in  which  case  the  last  N  —  M  rows  and 
columns  of  A  arc  assumed  to  be  zero.  In  the  QP  case,  H  should  not  be  greater  than  N; 
if  it  is,  the  last  M  —  N  rows  of  A  are  ignored. 

If  the  problem  is  specified  as  type  LSI,  LS2,  LS3  or  LS4,  M  is  also  the  dimension  of  the 
array  B.  Note  that  all  possibilities  (M  <  N,  M  =  N  and  H  >  N)  are  allowed. 

N  (Input)  The  number  of  variables,  i.e.,  the  dimension  of  X.  (N  nmst  be  positive.) 

NCLIN  (Input)  The  number  of  general  linear  constraints  in  the  problem,  (NCLIN  may  be 

zero.) 

NROWC  (Input)  The  declared  row  dimension  of  C.  (NROWC  must  be  at  least  1  and  at  least 

NCLIN.) 

NROWA  (Input)  The  declared  row  dimension  of  the  array  A.  (NROWA  must  be  at  least  1  and 

at  least  M.) 

C  (Input)  A  real  array  of  declared  dimension  (NROWC,*),  where  the  second  dimension 

must  be  at  le.'i.st  N.  The  i-th  row  of  C  contains  the  coefficients  of  the  i-th  general 
constraint,  t  =  1  to  NCLIN.  If  NCLIN  is  zero,  C  is  not  accessed;  the  actual  parameter 
may  then  be  any  convenient  array  or  an  array  with  dimension  (1,1). 

BL  (Input)  A  real  array  of  dimension  at  least  N  +  NCLIN  that  contains  the  lower  bounds 

for  all  the  constraints,  in  the  following  order  (which  is  also  observeti  for  BU,  ISTATE, 
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anil  CLAMDA);  the  lirst  M  I'lcinonts  of  BL  contain  the  lower  boniuls  on  the  variables;  if 
..'CLIN  ;•  0.  the  next  NCLIN  elements  of  BL  contain  the  lower  boniuls  for  the  general 
linear  eoustraints.  In  order  for  tlie  (iroblem  siiecilication  to  be  nieaniiigfiil,  it  is 
reijiiired  that  BL(j)  <  BU(7)  for  all  j.  To  specify  a  non-existent  lower  bound  (i.e., 
tj  '  -  X.),  the  value  used  must  satisfy  BL(j)  <  — BIGBND,  where  BIGBND  is  the  value  of 
the  o])tioiial  jiarameter  Infinite  Bound,  whose  default  value  is  lO'”  (see  Section  4.2). 
To  sju-eify  the  j-th  constr.aint  as  an  equality,  the  user  must  set  BL{j)  =  BU(j)  =  fi, 
say,  where  \ii\  <  BIGBND. 

(Input)  A  real  array  of  dimension  at  least  N  -t- NCLIN  that  contains  the  upper  bounds 
for  all  the  constraints,  in  the  same  order  descrilu'd  above  iinrler  BL.  To  sjiecify  a 
non-existent  ujiper  bound  (i.e.,  Uj  =  oo),  the  value  used  must  satisfy  BU(j/)  >  BIGBND. 

(Input)  A  real  array  of  dimension  at  least  N  containing  the  coefficients  of  the  explicit 
linear  term  of  the  objective  function.  If  the  problem  is  of  type  FP,  QPl,  QP3,  LSI  or 
LS3.  CVEC  is  not  accessed;  CVEC  may  then  be  declared  to  be  of  dimension  (1),  or  the 
actual  parameter  may  be  any  convenient  array. 

(Input)  An  integer  array  of  ilimension  at  least  N  +  NCLIN.  ISTATE  need  not  be 
iuitiali/.ed  if  Cold  Start  (the  default)  is  specified.  For  a  Warm  Start,  ISTATE  specifies 
the  desired  status  of  the  constraints  at  the  start  of  the  feasibility  phase.  Tlie  ordering 
of  ISTATE  is  the  same  as  that  described  above  for  BL,  i.e.,  the  first  N  components  of 
ISTATE  refer  to  the  upper  and  lower  bounds  on  the  variables,  and  components  N  +  1 
througli  N  -t-  NCLIN  refer  to  the  upper  and  lower  bounds  on  Cx.  Possible  values  for 
ISTATE  are: 

ISTATE(j)  Meaning 

0  The  corresponding  constr.aint  should  not  be  in  the  initial  working  set. 

1  Tlie  constraint  should  be  in  the  initial  working  set  at  its  lower  bound. 

2  The  constraint  should  be  in  the  initial  working  set  at  its  upper  bound. 

.3  Tlie  constraint  should  be  in  the  initial  working  set  as  an  equality.  This 

value  must  not  be  specified  unless  BL(j)  =  BU(j).  The  values  1,  2  or  3 
all  have  the  same  effect  when  BL(j)  =  BU(j). 

Other  values  of  ISTATE  are  also  acceptable.  In  particular,  if  LSSOL  has  been  called 
previously  with  the  same  values  of  N  and  NCLIN,  ISTATE  already  contains  satisfactory 
information. 

(Output)  If  LSSOL  exits  with  INFORM  =  0,  1  or  3,  the  values  in  the  array  ISTATE  in¬ 
dicate  the  status  of  the  constraints  in  the  active  set  at  the  solution.  Otherwise,  ISTATE 
indicati's  the  comimsition  of  the  working  set  at  the  final  iterate.  The  significance  of 
each  possible  value  of  ISTATE(j)  is  a.s  follows; 

ISTATE(j)  Meaning 

-2  The  constraint  violates  its  lower  bound  by  more  than  the  feasibility  tol¬ 
erance. 

-1  The  constraint  violates  its  upper  bound  by  more  than  the  feasibility 
tolerance. 

The  constraint  is  satisfied  to  within  the  feasibility  tolerance,  but  is  not 
in  the  working  set. 
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1  This  inequality  constraint  is  included  in  the  working  set  at  its  lower 
bound. 

2  This  inequality  constraint  is  included  in  the  working  set  at  its  upper 
bound. 

3  The  constraint  is  included  in  the  working  set  as  an  equality.  This  value 
of  ISTATE  can  occur  only  when  BL(j)  =  BU(j). 

KX  (Input)  An  integer  array  of  dimension  at  least  N.  KX  must  be  defined  on  input  for 

problems  QP3,  QP4,  LS3  or  LS4,  i.e.,  problems  in  which  A  is  specified  as  an  upper- 
trapezoidal  matrix.  KX  must  define  the  order  of  the  columns  of  the  matrix  A  with 
respect  to  the  ordering  of  X.  Thus,  if  KX(1)  =  5,  column  1  of  A  is  the  column  associated 
with  variable  X(5).  For  problems  of  type  FP,  LP,  QPl,  QP2,  LSI  or  LS2,  KX  need  not 
be  initialized. 

(Output)  KX  gives  the  order  of  the  columns  of  A  with  respect  to  the  ordering  of  X, 
as  described  above. 

X  (Input)  A  real  array  of  dimension  at  least  N.  X  contains  the  initial  estimate  of  the 

solution. 

(Output)  X  is  the  last  iterate  of  LSSOL.  If  INFORM  =  0,  1  or  3,  X  will  be  an  estimate 
of  the  solution. 

A  (Input)  A  real  array  of  dimension  (NROWA,*),  where  the  second  dimension  must  be 

at  least  N.  A  defines  the  data  matrix  A  in  LCLS. 

If  the  problem  is  of  type  FP  or  LP,  A  is  not  accessed  and  may  be  dimensioned  (1,1). 

If  the  problem  is  of  type  QPl  or  QP2,  the  first  M  rows  and  columns  of  A  must  contain 
the  leading  M  by  M  rows  and  columns  of  the  symmetric  Hessian  matrix.  Only  the 
diagonal  and  upper-triangular  elements  of  the  leading  M  rows  and  columns  of  A  ar^ 
referenced.  The  remaining  elements  arc  assumed  to  be  zero  and  need  not  be  assigned. 

For  problems  QP3,  QP4,  LS3  or  LS4,  the  first  M  rows  of  A  must  contain  an  M  by  N  upper- 
trapezoidal  factor  of  either  the  Hessian  matrix  or  the  least-squares  matrix,  ordered 
according  to  the  KX  array  (see  above).  The  factor  need  not  be  of  full  rank,  i.e.,  some  of 
the  diagonals  may  be  zero.  However,  as  a  general  rule,  the  larger  the  dimension  of  the 
leading  non-singular  submatrix  of  A,  the  fewer  iterations  will  be  required.  Elements 
outside  the  upper-triangular  part  of  the  first  M  rows  of  A  are  assumed  to  be  zero  and 
need  not  be  assigned. 

If  a  constrained  least-squares  problem  contains  a  very  large  number  of  observations, 
storage  limitations  may  prevent  storage  of  the  entire  least-squares  matrix.  In  such 
cases,  the  user  should  transform  the  original  A  into  a  triangular  matrix  before  the 
call  to  LSSOL  and  solve  the  problem  as  type  LS3  or  LS4. 

(Output)  If  the  problem  is  of  type  LS  or  QP,  A  contains  the  upper-triangular  matrix 
iZ  of  (7),  with  columns  ordered  as  indicated  by  KX  (see  above).  This  matrix  may 
be  used  to  obtain  the  variance-covariance  matrix  or  to  recover  the  u[)pcr-triangular 
factor  of  the  original  least-squares  matrix. 

B  (Input)  A  real  array  of  dimension  at  least  M.  If  the  problem  is  of  type  FP,  LP  or  QP, 

B  is  not  accessed  and  may  be  dimensioned  (1).  If  the  problem  is  of  type  LS,  B  must 
contain  the  vector  of  observations  b  in  problem  LCLS. 
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(Output)  On  exit  from  a  problem  of  type  LS,  B  contains  the  transformed  rcsidtial 
vector  (9). 

INFORM  (Output)  An  integer  that  indicates  the  result  of  LSSOL.  (If  Print  Level  >  0,  a 
short  description  of  INFORM  is  printed.)  The  possible  values  of  INFORM  are: 

INFORM  Meaning 

0  X  is  a  strong  local  minimum.  (The  projected  gradient  is  negligible,  the 
Lagrange  imdtipliers  are  optimal,  and  11^  is  non-singular.) 

1  X  is  a  weak  local  minimum.  (The  projected  gradient  is  negligible,  the 
Lagrange  multipliers  are  optimal,  but  is  singular  or  there  is  a  small 
multiplier.)  This  means  that  the  final  X  i.s  not  unique. 

2  The  solution  appears  to  be  unbounded.  This  value  of  INFORM  implies 
that  a  step  as  large  as  Infinite  Bound  would  have  to  be  taken  in  order 
to  continue  the  algorithm.  This  situation  can  occur  only  when  A  is 
singular,  there  is  an  explicit  bnear  term,  and  at  least  one  variable  has 
no  upper  or  lower  bound. 

3  No  feasible  point  was  found,  i.e.,  it  was  not  possible  to  satisfy  all  the 
constraints  to  within  the  fea.sibility  tolerance.  In  tliis  case,  the  constraint 
violations  at  the  final  X  will  reveal  a  value  of  the  tolerance  for  which  a 
feasible  point  will  exist — for  example,  if  the  feasibility  tolerance  for  each 
violated  constraint  exceeds  its  Residual  at  the  final  point.  The  modified 
problem  (with  an  altered  feasibility  tolerance)  may  then  be  solved  using 
a  Warm  Start  (see  Section  4). 

4  The  limiting  number  of  iterations  (determined  by  the  parameters  Feasi¬ 
bility  Phase  Iterations  and  Optimality  Phase  Iterations)  was 
reached  before  normal  termination  occurred. 

5  The  algorithm  could  be  cycUng,  since  a  total  of  50  changes  were  made 
to  the  working  set  without  altering  X. 

6  An  input  parameter  is  invalid. 

ITER  (Output)  An  integer  that  gives  the  total  number  of  iterations  performed  in  the 

feasibility  phase  and  the  optimality  phase. 

OBJ  (Output)  The  value  of  the  objective  function  at  X  if  X  is  feasible,  or  the  sum  of 

infeasibilities  at  X  otherwise.  If  the  problem  is  of  type  FP  and  X  is  feasible,  OBJ  is  zero. 

CLAMDA  (Output)  A  real  array  of  dimension  at  least  N  -f  NCLIN  that  contains  the  Lagrange 
multiplier  for  every  constraint  with  respect  to  the  current  working  set.  The  ordering 
of  CLAMDA  follows  the  convention  given  above  under  BL,  i.e.,  the  first  N  components 
contain  the  multipliers  for  the  bound  constraints  on  the  variables,  and  the  remaining 
components  contain  the  multipliers  for  the  general  linear  constraints.  If  ISTATE(_;)  =  0 
(i.e.,  constraint  j  is  not  in  the  working  set),  CLAMDA(j)  is  zero.  If  X  is  optimal, 
CLAHDA(j)  should  be  non-negative  if  ISTATE(j)  =  1  and  non-positive  if  ISTATE(j)  =  2. 

3.2.  Workspace  parameters 

IW  (Input)  An  integer  array  of  dimension  LENIW  that  provides  integer  workspace  for 

LSSOL. 
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LENIW  (Input)  The  dimension  of  IW.  LENIW  must  be  at  least  N. 

W  (Input)  A  real  array  of  dimension  LENW  that  provides  real  workspace  for  LSSOL. 

LENW  (Input)  The  dimension  of  W.  If  the  problem  is  of  type  FP  and  N  <  NCLIN,  LENW  must 

be  at  least  2  N^  +  6  N  +  6  NCLIN.  If  the  problem  is  of  type  FP  and  0  <  NCLIN  <  N,  LENW 
must  be  at  least  2  (NCLIN  +  1)*  +  6N  +  6NCLIN.  If  NCLIN  =  0,  LENW  must  be  at  least 
6N. 

If  the  problem  is  of  type  LP  and  N  <  NCLIN,  LENW  must  be  at  least  2N^  +7N  -f- 6 NCLIN. 
If  the  problem  is  of  type  LP  and  N  >  NCLIN  >  0,  LENW  must  be  at  least  2  (NCLIN  + 
1)^  +  7  N  +  6  NCLIN.  If  the  problem  is  of  type  LP  and  NCLIN  =  0,  LENW  must  be  at  least 
7N. 

For  problems  QPl,  qP3,  LSI  and  LS3,  LENW  must  be  at  least  2N^  +  9N  +  6NCLIN  if 
NCLIN  >  0,  and  at  least  9  N  if  NCLIN  =  0.  For  problems  qP2,  qP4,  LS2  and  LS4,  LENW 
must  be  at  least  2N^  +  ION  4-  6  NCLIN  if  NCLIN  >  0,  and  at  least  ION  if  NCLIN  =  0. 

If  Print  Level  >  0,  the  amounts  of  workspace  provided  and  required  are  printed.  As  an  alterna¬ 
tive  to  computing  LENIW  and  LENW  from  the  formulas  given  above,  the  user  may  prefer  to  obtain 
appropriate  values  from  the  output  of  a  preliminary  run  with  a  positive  value  of  Print  Level  and 
LENIW  and  LENW  set  to  1.  (LSSOL  will  then  terminate  with  INFORM  =  6.) 
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4.  OPTIONAL  INPUT  PARAMETERS 

Several  optional  parameters  in  LSSOL  define  choices  in  the  problem  specification  or  the  algorithm 
logic.  In  order  to  reduce  the  number  of  formal  parameters  of  LSSOL,  these  optional  parameters 
have  a.ssociated  default  values  (see  Section  4.2)  that  arc  appropriate  for  most  problems.  Therefore, 
the  user  need  specify  only  those  parameters  whose  values  are  to  be  different  from  their  default 
values.  The  remainder  of  this  section  can  be  skippeil  by  users  who  wish  to  use  the  default  values 
for  all  optional  parameters. 

Each  optional  parameter  is  defined  by  a  single  character  string  of  up  to  72  characters,  con¬ 
taining  one  or  more  items.  The  items  associated  with  a  given  option  must  be  separated  by  spaces 
or  equal  signs  (=).  Alphabetic  characters  may  be  upper  or  lower  case.  An  example  of  an  optional 
parameter  is  the  string 

Print  level  =  5 

For  each  option,  the  string  contains  the  following  items. 

1.  The  keyword  (required  for  all  options). 

2.  A  phrase  (one  or  two  words)  that  qualifies  the  keyword  (only  for  some  options). 

3.  A  riumher  that  specifies  either  an  INTEGER  or  a  REAL  value  (only  for  some  options). 

Such  numbers  may  be  up  to  16  contiguous  characters  in  Fortran  77’s  I,  F,  E  or  D 
formats,  terminated  by  a  space. 

Blank  strings  and  comments  are  ignored  and  may  be  used  to  improve  readability.  A  coinmcnt  begins 
with  an  asterisk  (♦)  and  all  subsequent  characters  are  ignored.  If  the  string  is  not  a  comment  and 
is  not  recognized,  a  warning  message  is  printed  on  the  specified  output  device  (see  Section  7.5). 
Synonyms  are  recognized  for  some  of  the  keywords,  and  abbreviations  may  be  used. 

The  following  arc  examples  of  valid  option  strings  for  LSSOL; 

NOLIST 
warm  start 
COLD  START 

Problem  type  =  Least  Squares 
Problem  type  =  LP 
Problem  Type  QP4 

Feasibility  toleramce  l.OE-8  *  for  IBM  in  double  precision 
CRASH  TOLERANCE  =  .002 

•  This  string  will  be  completely  ignored. 

Feasibility  phase  iteration  limit  100 
Optimality  phase  iteration  limit  =  10  ♦ 

4.1.  Specification  of  the  optional  parameters 
Optional  parameters  may  be  specified  in  two  ways,  as  follows. 

•  Using  subroutine  LSFILE  and  an  external  file 

The  subroutine  LSFILE  provided  with  the  LSSOL  package  will  read  options  from  an  external  options 
file,  and  shoidd  be  called  before  a  call  to  LSSOL.  Each  line  of  the  options  file  defines  a  single  optional 
parameter.  The  file  must  begin  with  Begin  and  end  with  End.  (An  options  file  consisting  only  of 
these  two  lines  corresponds  to  supplying  no  options.) 

The  specification  of  LSFILE  is 

SUBROUTINE  LSFILEf  lOPTNS,  INFORM  ) 

INTEGER  lOPTNS,  INFORM 


OPTIONAL  INPUT  PARAMETERS 


IS 


IQPTNS  must  be  the  unit  number  of  the  options  file,  in  the  riinge  (0,99),  and  is  unchanged  on  exit 
from  LSFILE.  INFORM  need  not  be  set  on  entry.  On  '•eturn.  INFORM  will  be  0  if  the  file  is  a  valid 
options  file  ami  lOPTNS  is  in  the  correct  rang»'.  INFORM  will  be  set  to  1  if  lOPTNS  is  out  of  range, 
and  will  be  set  to  2  if  the  file  does  not  begin  with  Begin  or  end  with  End. 

An  example  of  a  valid  options  file  is 

Begin 

Print  level  =  B 
Problem  type  LP 
End 

If  the  options  file  is  on  unit  number  5,  it  can  be  read  by  the  call 

CALL  LSFILEC  5.  INFORM  ) 


•  Using  subroutine  LSOPTN 

The  second  method  of  setting  the  optional  parameters  is  through  a  series  of  calls  to  the  subroutine 
LSOPTN  provided  with  the  LSSOL  package.  The  specification  of  LSOPTN  is 

SUBROUTINE  LS0PTN(  STRING  ) 

CHARACTER* (*)  STRING 

STRING  must  be  a  single  valid  option  string  (see  above),  and  will  be  unchanged  on  exit.  LSOPTN 
must  be  called  once  for  every  optional  parameter  to  be  set.  An  example  of  a  call  to  LSOPTN  is 

CALL  LSOPTN (  'Print  level  =  6'  ) 


•  Use  of  the  Nolist  and  Defaults  option 

In  general,  each  user-specified  optional  parameter  is  printed  as  it  is  read  or  defined.  By  using  the 
special  parameter  Nolist,  the  user  may  suppress  this  printing  for  a  given  call  of  LSSOL.  To  take 
effect,  Nolist  must  be  the  first  parameter  specified  in  the  options  file;  for  example. 

Begin 

Nolist 

Problem  type  LP 
End 

Alternatively,  the  first  call  to  LSOPTN,  before  or  after  a  call  to  LSSOL,  must  be 

CALL  LSOPTNC  'Nolist'  ). 

All  parameters  not  specified  by  the  user  are  automatically  set  to  their  default  values.  Any 
optional  parameters  that  are  set  by  the  user  are  not  altered  by  LSSOL,  and  hence  changes  to  the 
options  are  cumii/ative.  For  example,  calling  LS0PTN(  'Print  level  =  5’  )  sets  the  print  level 
to  5  for  all  subsequent  calls  to  LSSOL  until  it  is  reset  by  the  user.  The  only  exception  to  this 
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nil*'  is  pcrinittcd  by  the  special  optional  parameter  Defaults,  whoso  effect  is  to  reset  nil  optional 
parameters  to  their  default  values.  For  example,  in  the  following  situation 

CALL  LSSOL  (  . . .  ) 

C 

CALL  LSOPTN(  'Print  level  5'  ) 

CALL  LSQPTNC  'Iteration  limit  =  100'  ) 

CALL  LSSOL  (  . . .  ) 

C 

CALL  LSOPTNC  'Defaults'  ) 

CALL  LSSOL  (  . . .  ) 

the  first  and  last  runs  of  LSSOL  will  occur  with  the  default  parameter  settings,  but  in  the  second 
run,  the  print  level  and  iteration  limit  are  altered. 

4.2.  Description  of  the  optional  parameters 

The  following  list  (in  alphabetical  order)  gives  the  valid  options.  For  each  option,  we  give  the 
keyword,  any  essential  ojitional  qualifiers,  the  <lefault  value,  and  the  definition.  The  miniiimin 
abbreviation  of  each  keyword  is  underlined.  If  no  characters  of  an  optional  qualifier  are  uiulerlined, 
the  qualifer  may  be  omitted.  The  letter  a  denotes  a  phrase  (character  string)  that  cjiialifies  an 
option.  The  letters  i  and  r  denote  INTEGER  and  REAL  values  required  with  certain  options.  The 
number  c  is  a  generic  notation  for  machine  precision. 

Cold  Start  Default  =  Cold  Start 

Warm  Start 

This  oi)tion  specifies  how  the  initial  working  set  is  chosen.  With  a  cold  start,  LSSOL  chooses 
the  initial  working  set  based  on  the  values  of  the  variables  and  constraints  at  the  initial  point. 
Broadly  speaking,  the  initial  working  set  will  include  equality  constraints  and  bounds  or  inequality 
constraints  that  violate  or  “nearly’’  satisfy  their  bounds  (to  within  Crash  Tolerance;  see  below). 

With  a  warm  start,  the  user  must  provide  a  valid  definition  t)f  every  clement  of  the  array 
ISTATE  (see  Section  3  for  the  definition  of  this  array).  LSSOL  will  override  the  user's  specification 
of  ISTATE  if  necessary,  so  that  a  poor  choice  of  the  working  set  will  not  cause  a  fatal  error.  A  warm 
start  will  be  advantageous  if  a  good  estimate  of  the  initial  working  set  is  available  -for  exantple, 
when  LSSOL  is  called  repeatedly  to  solve  related  problems. 

Crash  Tolerance  r  Default  =  .01 

This  value  is  used  in  conjunction  with  the  optional  parameter  Cold  Start  (the  default  value)  when 
LSSOL  selects  an  initial  working  set.  If  0  <  r  <  1,  the  initial  working  set  will  inclmle  bounds  or 
general  inequality  constraints  that  lie  within  r  of  their  bounds.  In  particidar,  a  constraint  of  the 
form  cj'x  >  I  will  be  included  in  the  initial  working  set  if  |cji  —  <  r(l-f|f|).  Ifr  <  0  or  r  >  1, 

the  defa)dt  valtie  is  used. 

Feasibility  Phase  Iteration  Limit  Default  =  inax(50,  5(n mt )) 

Optimality  Phase  Iteration  Limit  »2  Dcfaidt  =  max(50, 5(n -f- mt )) 

The  scalars  f'l  and  *2  specify  the  maximum  number  of  iterations  allowed  in  the  feasibility  and  opti¬ 
mality  phases.  Optimality  Phase  Iteration  Limit  is  equivalent  to  Iteration  Limit.  Setting 
t'l  =  0  and  Print  Level  >  0  means  that  the  workspace  needed  will  be  computed  and  printed,  but 
no  iterations  will  be  performed. 
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Feasibility  Toleremce  r  Default  =  i/c 

If  r  >  0,  r  drfiiu's  tlie  luaxiniiini  acceptable  nhsoliitc  violati«)n  in  each  constraint  at  a  “feasible” 
point;  i.e.,  a  constraint  is  considered  satisfled  if  its  violation  docs  not  exceed  r.  For  example,  if  the 
variables  and  the  coefficients  in  the  general  constraints  are  of  order  unity,  and  the  latter  arc  correct 
to  about  6  decimal  digits,  it  would  be  appropriate  to  specify  r  as  10~®.  If  r  <  0,  the  default  value 
is  used. 


Infinite  Bound  Size 


Default  = 


If  r  >  0.  r  defines  the  “infinite”  bound  BIGBND  in  the  definition  of  the  problem  constraints.  Any 
upper  bound  greater  than  or  equal  to  BIGBND  will  be  regarded  as  plus  infinity  (and  similarly  for  a 
U)wer  bound  less  than  or  equal  to  —BIGBND).  If  r  <  0,  the  default  value  is  used. 

Infinite  Step  Size  r  Default  =  max(BIGBND,  10^®) 

If  r  >  0,  r  specifies  the  magnitude  of  the  change  in  variables  that  will  be  considered  a  step  to 
an  tinbounded  solution.  (Note  that  an  unbounded  solution  can  occur  only  when  the  Hessian  is 
singular  and  the  objective  contains  an  explicit  linear  term.)  If  the  change  in  x  during  an  iteration 
would  exceed  the  value  of  Infinite  Step,  the  objective  function  is  considered  to  be  unbounded 
below  in  the  feasible  region.  If  r  <  0,  the  default  value  is  used. 


Iteration  Limit  i 

Iters 

Itns 

See  Feasibility  Phase  Iteration  Limit  above. 

Optimality  Phase  Iteration  Limit  i 

See  Feasibility  Phase  Iteration  Limit  above. 

Print  Level  » 


Default  =  max(50, 5(n  +  nit,)) 


Default  =  max(50, 5(n  +  rui,)) 


Default  =  10 


The  value  of  i  controls  the  amount  of  printout  produced  by  LSSOL,  as  indicated  below. 

t  Output 

0  No  output. 

1  The  final  solution  only. 

5  One  line  of  output  for  each  iteration  (no  printout  of  the  final  solution). 

>  10  The  final  solution  and  one  line  of  output  for  each  iteration. 

>  20  At  each  iteration,  the  Lagrange  multipliers,  the  variables  x,  the  constraint 

values  Cx  and  the  constraint  status. 

>  30  At  each  iteration,  the  diagonal  elements  of  the  matrix  T  associated  with  the  TQ 

factorization  (3)  of  the  working  set,  and  the  diagonal  elements  of  the  triangular 
matrix  R. 

Problem  Type  o  Default  =  LSI 

This  option  specifies  the  type  of  o  ^ active  function  to  be  minimized  dtiring  the  optimality  phase. 
The  following  arc  the  ten  optional  keywords  and  the  dimensions  of  the  arrays  that  imist  be  specified 


■.v.v, 
•  V  V 
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to  define  the  objective  function: 

FP  A,  B  and  CVEC  not  accessed; 

LP  A  and  B  not  accessed,  CVEC(M); 

QPl  A(NROWA,N)  symmetric,  B  and  CVEC  not  referenced; 

QP2  A(NROWA,N)  symmetric,  B  not  referenced,  CVEC(N); 

QP3  A(NROWA,N)  upper-trapezoidal,  KX(N),  B  and  CVEC  not  referenced; 

QP4  A(NR0WA,N)  uj)per-trapezoidal,  KX(N),  B  not  referenced,  CVEC(N); 

LSI  A ( NROWA ,  N ) ,  B  ( M ) ,  CVEC  not  referenced; 

LS2  A(NROWA,W),  B(M),  CVEC(N); 

LS3  A(NROWA,N)  upper-trapezoidal,  KX(N),  B(M),  CVEC  not  referenced; 

LS4  A(NRQWA,N)  upper-trapezoidal,  KX(N),  B(M),  CVEC(N). 

The  options  Least  Squares  and  LSQ  are  equivalent  to  the  default  option  LSI.  The  options 
Linear  program  and  Quadratic  program  are  equivalent  to  LP  and  QP2  respectively.  If  =  0,  i.e., 
the  objective  function  is  purely  linear,  the  efficiency  of  LSSOL  may  be  increased  by  specifying  a  as 
LP  (or  Linear  Program). 

R^k  Tolerance  r  Default  =  y/e 

If  0  <  r  <  1,  r  enables  the  user  to  control  the  estimation  of  the  rank  of  A  and  the  triangular  factor 
i?i  (see  Section  2).  If  pi  denotes  the  function  pi  =  max{|7Zii |,  |/?22|>  •  •  •  i  the  rank  of  R  is 

defined  to  be  smallest  index  i  such  that  <  T\pi^i\.  If  r  <  0  or  r  >  1,  the  default  value  is 

used. 


4.3.  Optional  parameter  checklist  and  default  values 

For  easy  reference,  the  following  sample  LSOPTN  Ust  shows  all  vaUd  keywords  and  their  default 
values.  The  default  options  Feasibility  Tolerance  and  Rank  Tolerance  depend  upon  e,  the 
relative  precision  of  the  machine  being  used.  The  values  given  here  correspond  to  double  precision 
arithmetic  on  IBM  360  and  370  systems  and  their  successors  (e  %  2.22  x  10  ^®).  Similar  values 
would  apply  to  any  machine  having  about  16  decimal  digits  of  precision. 

*  List  of  optional  parameters. 

Cold  Start 
Crash  Tolerance 
Feasibility  Toleremce 
Infinite  Bound 
Infinite  Step 

Feasibility  Phase  Iteration  Limit 
Optimality  Phase  Iteration  Limit 
Print  Level 
Problem  Type 
Rank  Tolerance 


.01 

l.lE-8 

l.OE+10 

l.OE+10 

50 

50 

10 

Least  squares 
l.lE-8 


* 

* 

♦  Plus  infinity 

♦  or  5(n  -1-  thl) 

*  or  5(n  -f-  ttil) 

* 

*  or  LSI 

♦  y/e 


5.  DESCRIPTION  OF  THE  PRINTED  OUTPUT 
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5.  DESCRIPTION  OF  THE  PRINTED  OUTPUT 

This  section  describes  the  intermediate  printout  produced  by  LSSOL.  To  aid  interpretation  of  the 
printed  results,  we  repeat  the  convention  for  nuiid>ering  the  constraints:  indices  1  through  N  refer  to 
the  bounds  on  the  variables,  and  indices  N  +  1  through  N  +  NCLIN  refer  to  the  general  constrsiints. 
When  the  status  of  a  constraint  changes,  the  index  of  the  constraint  is  printed,  along  with  the 
designation  “L”  (lower  bound),  “U”  (upper  bound),  “E”  (equality),  “T”  (temporary  bound)  or  “Z” 
(artificial  constraint). 

When  Print  Level  >  5,  the  following  line  of  output  is  produced  at  every  iteration.  In  all 
cases,  the  values  of  the  quantities  printed  are  those  in  effect  on  completion  of  the  given  iteration. 
Itn  is  the  iteration  count. 

Jdel  is  the  index  of  the  constraint  deleted  from  the  working  set.  If  Jdel  is  zero,  no 

constraint  was  deleted. 

Jadd  is  the  index  of  the  constraint  added  to  the  working  set.  If  Jadd  is  zero,  no 

constraint  was  added. 

Stop  is  the  step  taken  along  the  computed  search  direction.  If  a  constraint  is  added 

during  the  current  iteration  (i.e.,  Jadd  is  positive),  Step  will  be  the  step  to  the 
nearest  constraint.  During  the  optimabty  phase,  the  step  can  be  greater  than 
one  only  if  the  factor  is  singular. 

Minf  is  the  number  of  violated  constraints  (infeasibilities).  This  number  will  be  zero 

during  the  optimality  phase. 

Sinf /Objective  is  the  value  of  the  current  objective  function.  If  X  is  not  feasible,  Sinf  gives 
a  weighted  sum  of  the  magnitudes  of  constraint  violations.  If  X  is  feasible, 
Objective  is  the  value  of  the  objective  function  of  LCLS.  The  output  line  for 
the  final  iteration  of  the  feasibility  phase  (i.e.,  the  first  iteration  for  which  NINE 
is  zero)  will  give  the  viduc  of  the  true  objective  at  the  first  feasible  point. 
During  the  optimality  phase,  the  value  of  the  objective  function  will  be  non- 
increasing.  During  the  feasibility  phase,  the  number  of  constraint  infeasibilities 
will  not  increase  until  either  a  feasible  point  is  found,  or  the  optimality  of  the 
multipliers  implies  that  no  feasible  point  exists.  Once  optimal  multipliers  are 
obtained,  the  number  of  infcasibilities  can  increase,  but  the  sum  of  infeasi¬ 
bilities  will  cither  remain  constant  or  be  reduced  until  the  minimum  sum  of 
infcasibilities  is  found. 

End  is  the  number  of  simple  bound  constraints  in  the  current  working  set. 

Lin  is  the  number  of  general  linear  constraints  in  the  current  working  set. 

Nz  is  the  number  of  columns  of  Z  (see  Section  2),  The  value  of  Nz  is  the  number 

of  variables  minus  the  number  of  constraints  in  the  working  set;  i.e.,  Mz  = 
N-  (End -t- Lin).  A  zero  value  of  Nz  implies  that  x  lies  at  a  vertex  of  the  feasible 
region. 

Nzl  is  the  number  of  columns  of  Zi  (see  Section  2).  Nzl  is  the  dimension  of  the 

subspace  in  which  the  objective  function  is  currently  being  minimized.  If  Nzl 
is  less  than  Nz,  the  current  is  singular. 

Norm  Gf  is  the  Euclidean  norm  of  the  gradient  of  the  objective  function  with  respect  to 

the  free  variables,  i.e.,  variables  not  currently  held  at  a  bound. 

Norm  Gzl  is  the  Euclidean  norm  of  the  projected  gradient  with  respect  to  Z\. 

During  the  optimality  phase,  this  norm  will  be  approximately  zero  after  a  unit 
step. 
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Cond  T 
Cond  Rzl 


When  Print 
of  LSSOL  includes 
are  assigned  to  all 
The  following 
Variable 
State 


Value 

Lower  bound 


Upper  bound 


Lagr  multiplier 


Residual 


is  a  lower  bound  on  the  condition  number  of  the  working  set. 
is  a  lower  bound  on  the  condition  number  of  the  triangular  factor  R-^  (the  first 
Nzl  rows  and  columns  of  the  factor  R^.  If  the  problem  is  specified  to  be  of 
type  LP,  or  the  estimated  rank  of  the  data  matrix  A  is  zero,  Cond  Rzl  is  not 
printed. 

Level  =  1  or  Print  Level  >  10,  the  summary  printout  at  the  end  of  execution 
a  listing  of  the  status  of  every  variable  and  constraint.  Note  that  default  names 
variables  and  constraints, 
describes  the  printout  for  each  variable. 

gives  the  name  (VARBL)  and  index  j  (j  =  1  to  N)  of  the  variable. 

gives  the  state  of  the  variable  (FR  if  neither  bound  is  in  the  working  set,  EQ  if 
a  fixed  variable,  LL  if  on  its  lower  bound,  UL  if  on  its  upper  bound).  If  Value 
lies  outside  the  upper  or  lower  bounds  by  more  than  the  feasibility  tolerance, 
State  will  be  “++”  or  “ — ”  respectively. 

is  the  value  of  the  variable  at  the  final  iteration. 

is  the  lower  bound  specified  for  the  variable.  (“None”  indicates  that  BL(ji)  < 
-BIGBND.) 

is  the  upper  bound  specified  for  the  variable.  (“None”  indicates  that  BU(7)  > 
BIGBND.) 

is  the  value  of  the  Lagrange  multiplier  for  the  associated  bound  constraint.  This 
will  be  zero  if  State  is  FR.  If  X  is  optimal,  the  multiplier  should  be  non-negative 
if  State  is  LL,  and  non-positive  if  State  is  UL. 

is  the  difference  between  the  variable  “Value”  and  the  nearer  of  its  bounds 
BL(j)  and  BU(j). 


The  meaning  of  the  printout  for  general  constraints  is  the  same  as  that  given  above  for  vari¬ 
ables,  with  “variable”  replaced  by  “constraint”,  with  the  following  change  in  the  heading; 

Linear  constr  is  the  name  (LNCON)  and  index  i  (i  =  1  to  NCLIN)  of  the  constraint. 
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6.  ERROR  RECOVERY 

Termination  Recommended  Action 

A  siiif'lc  underflow  will  always  occur  if  inacliine  constants  are  computed  automat¬ 
ically  (as  in  the  distributed  version  of  LSSOL;  sec  Section  7).  Other  floating-point 
underflows  may  occur  occasionally,  but  can  usually  be  ignored. 

If  the  printed  (uitput  before  the  overflow  error  contains  a  warning  about  serious 
ill-conditioning  in  the  working  set  when  adding  the  ji-th  constraint,  it  may  be  pos¬ 
sible  to  avoid  the  difKculty  by  increasing  the  magnitude  of  the  optional  i)arameter 
Feasibility  Tolerance  and  rerunning  the  program.  If  the  message  recurs  even 
after  this  change,  the  offending  linearly  dependent  constraint  (with  index  “y”) 
must  be  removed  from  the  problem.  If  a  warning  message  did  not  precede  the 
fatal  overflow,  contact  the  authors  at  Stanford  University. 

LSSOL  has  terminated  without  finding  a  feasible  point,  which  means  that  no  fea¬ 
sible  point  exists  for  the  given  feasibility  tolerance.  The  user  should  check  that 
there  are  no  constraint  redundancies.  If  the  data  for  the  constraints  are  accurate 
only  to  the  absolute  preci.sion  cr,  the  user  should  ensure  that  the  value  of  the  op¬ 
tional  parameter  Feasibility  Tolerance  is  greater  than  a.  For  example,  if  all 
elements  of  C  are  of  order  unity  and  are  accurate  only  to  three  <lecimal  places,  the 
optional  parameter  Feasibility  Tolerance  should  be  at  least  10 

The  value  of  the  optional  parameter  Iteration  Limit  may  be  too  small.  If  the 
method  appears  to  be  making  progress  (e.g.,  the  objective  function  is  being  sat¬ 
isfactorily  reduced),  increase  the  iterations  limit  and  rerun  LSSOL  (possibly  using 
the  warm  start  facility  to  specify  the  initial  working  set).  If  the  iteration  limit  is 
already  large,  btit  some  of  the  constraints  could  be  nearly  linearly  dependent,  check 
the  output  for  a  repeated  pattern  of  constraints  entering  and  leaving  the  working 
set.  (Near-dependencies  are  often  indicated  by  wide  variations  in  size  in  the  di¬ 
agonal  elements  of  the  T  matrix,  which  will  be  printed  if  Print  Level  >  30.)  In 
this  case,  the  algorithm  could  be  cycling  (see  the  comments  for  INFORM  =  5). 

INFORM  =  5  This  value  will  occur  if  50  iterations  are  performed  without  changing  X.  The  user 
should  check  the  printed  output  for  a  repeated  pattern  of  constraint  deletions  and 
additions.  If  a  sequence  of  constraint  changes  is  being  rei>eated,  the  iterates  are 
probably  cycling.  (LSSOL  does  not  contain  a  method  that  is  guaranteed  to  avoid 
cycling;  such  a  method  would  be  combinatorial  in  nature.)  Cycling  may  occur  in 
two  circumstances;  at  a  constrained  stationary  point  where  there  are  some  small 
or  zero  Lagrange  multipliers;  or  at  a  point  (usually  a  vertex)  where  the  constraints 
that  are  satisfled  exactly  are  nearly  linearly  dependent.  In  the  latter  case,  the  user 
has  the  option  of  identifying  the  offending  dependent  constraints  and  removing 
them  from  the  problem,  or  restarting  the  run  with  a  larger  value  of  the  optional 
parameter  Feasibility  Tolerance.  If  LSSOL  terminates  with  INFORM  -  5,  but 
no  suspicio\is  pattern  of  constraint  changes  can  be  observeil,  it  may  be  worthwhile 
to  restart  with  the  final  X  (with  or  without  the  warm  start  option). 


Underflow 

Overflow 


INFORM  =  3 


INFORM  =  4 
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7.  IMPLEMENTATIOIM  INFORMATION 

7.1.  Format  of  the  distribution  tape 

The  source  code  and  exain])le  proj^raui  for  LS.SOf.  are  distril)ut<‘d  on  a  magnetic  tape  containing  7 
files.  'I'lie  tape  characteristics  are  described  in  a  domiment  accompanying  the  tape;  normally  they 
are  9  track,  IGOO  bpi,  uidabeled,  ASCII,  80-character  reconls  {card  images),  dSOO-character  blocks. 

TIk'  following  is  a  list  of  the  files  and  a  summary  of  their  contents.  For  reference  purposes  wc 
give  a  name  to  each  file.  However,  the  names  will  not  be  recorded  on  unlabeled  tapes.  The  MACH 
and  LSCODE  files  are  composed  of  several  smaller  source  files  described  in  Section  7.3. 


i/e  Name 

Type 

Cardsf 

Description 

1.  DPMACH 

FORTRAN 

450 

Double- precision  source  file  1:  MCSUBS 

2.  DPLSCODE 

FORTRAN 

8250 

Double-i)recision  source  files  2  -5;  BLAS, 

3.  DPLSMAIW 

FORTRAN 

260 

Double-precision  source  file  LSMAIN 

4.  LSMAIN 

DATA 

6 

Options  file  for  LSMAIN 

5.  SPMACH 

FORTRAN 

450 

Single-precision  source  file  1 

6.  SPLSCODE 

FORTRAN 

8250 

Single-precision  source  files  2  5 

7.  SPLSHAIN 

FORTRAN 

260 

Single-precision  version  of  file  3 

t  Approximate  figure. 


One  MACH  and  one  LSCODE  file  shotdd  be  selected  for  any  given  installation.  DPMACH  and 
DPLSCODEare  inteiuled  for  machines  that  generally  require  doubh'  precision  com])utation.  Examples 
include  HIM  Systems  360,  370,  3033,  3081,  etc.;  Aitidahl  470,  Facom,  Fujitsu,  Hitachi,  and  other 
systems  analogous  to  IBM;  DEC  VAX  systems;  Data  General  MV/8000;  ICL  2900  series;  recent 
PRIME  systems;  DEC  Systems  10  and  20;  Honeywell  systems;  and  the  Univac  1100  scries. 

SPMACH  and  SPLSCODE  are  intended  for  machines  for  which  single  precision  is  suitably  accurate 
for  nui/ierical  computation.  Examples  include  the  Burrouglis  6700  and  7700  series;  the  CDC  6000 
and  7000  series  and  their  Cyber  counterparts;  and  the  Cray-1. 


7.2.  Installation  procedure 

1.  flbtain  the  appropriate  MACH  and  LSCODE  files  from  the  tape. 

2.  If  necessary,  edit  the  subroutine  MCHPAR  according  to  Section  7.5. 

3.  Decide  whether  or  not  to  split  the  LSCODE  file  into  files  BLAS  through  OPSUBS  as  suggested  in 
Section  7.3. 

4.  Compile  all  the  routines  that  were  originally  in  the  LSCODE  files  together  with  those  from  MACH. 
Run  them  in  conjunction  with  the  main  program  LSMAIN  from  either  File  3  or  File  7  and  the 
options  given  in  file  LSMAIN  DATA.  Check  the  output  against  that  shown  in  Section  8. 


7.3.  Source  files 

LSSOL  has  been  written  in  ANSI  (1977)  Fortran  and  tested  on  an  IBM  3n8lK  computer  using  the 
IBM  Fortran  77  compiler  V'S  Fortran.  Certain  unavoidable  machine  dependencies  are  confined  to 
the  routine  MCHPAR. 

The  source  code  is  divided  into  5  logical  parts.  For  ease  of  handling,  these  are  combined  into 
the  MACH  and  LSCODE  files  on  the  distribution  tape,  but  for  subsequent  maintenance  we  recommend 
that  5  sci)arate  files  be  kept.  In  the  description  below  we  suggest  a  name  for  eacli  file  and  summarize 
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its  purpose.  Wc  then  list  the  names  of  the  Fortran  subroutines  and  functions  involved.  The  naming 
convention  used  should  minimize  the  risk  of  a  clash  with  user-written  routines. 


File  1. 


MCSUBS  Computes  machine-dependent  constants. 
MCHP-iR  MCEPS  MCENVl  MCENV2  MCSTOR 


File  2.  BLAS  Basic  Linear  Algebra  Subprograms  (a  subset). 

DASUM  DAXPY  DCOPY  DDOT  DWRM2  DSWAP  DSCAL  IDAMAX 

These  routines  are  functionally  similar  to  members  of  the  BLAS  package  (Lawson  et  al., 
1979).  If  possible  they  should  be  replaced  by  authentic  BLAS  routines.  Versions  may 
exist  that  have  been  tuned  to  your  particular  machine. 

DGEMV  DGERl 

These  routines  are  functionally  similar  to  members  of  the  Level  2  BLAS  packages  (Don- 
garra  et  al,  1985). 

DCQND  DDIV  DDSCL  DLOAD  DNORM  DSSQ  DSWAP  ICOPY 

IDRANK  ILOAD 

These  are  additional  utility  routines  that  could  be  tuned  to  your  machine.  DLOAD  is  used 
the  most  frequently,  to  load  a  vector  with  a  constant  value. 

DR0T3  DR0T3G  DGEAPQ  DGEQR  DGEQRP  DGRFG 

These  linear  algebra  routines  arc  used  to  compute  and  update  various  matrix  factoriza¬ 
tions  in  LSSOL. 


File  3. 

CMSUBS 

General 

utility  routines. 

CMALF 

CMRIHD 

CMALFl 

CHTSOL 

CMCHK  CMFEAS 

CMPRT 

CMQMUL 

CMRSOL 

CMRSWP 

File  4. 

LSSUBS 

Least-squares  routines. 

LSADD 

LSFEAS 

LSOPTN 

LSADDS 

LSFILE 

LSPRT 

LSBNDS  LSCHOL 

LSGETP  LSGSET 

LSSETX  LSSOL 

LSCORE 

LSKEY 

LSCRSH 

LSLOC 

LSDEL 

LSMOVE 

LSDFLT 

LSHULS 

File  5. 

OPSUBS 

Option  string  handling  routines. 

OPFILE 

OPLOOK 

OPNUM  OPSCAN 

OPTOKN 

OPUPPR 

7.4.  Common  blocks 

Certain  Fortran  COMMON  blocks  are  used  in  the  LSSOL  source  code  to  communicate  between  sub¬ 
routines.  Their  names  are  listed  below. 

CMDEBG  LSDEBG  LSPARl  LSPAR2  SOLICH  S0L3CM  S0L4CM  S0L5CM 

S0L6CM  SOLMCH  SOLILS  S0L3LS 


7.5.  Machine-dependent  subroutines 

The  routine  MCHPAR  in  the  MACH  file  may  require  modification  to  suit  a  particular  machine  or  a 
non-standard  application. 


m 
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At  the  hi't'inuiii!;  of  LSSOL,  MCHPAR  is  called  to  assign  the  inachiue-dejxMwleiit  constants  and 
the  standard  inpAit  and  imlput  unit  innnbers.  These  i)aranieters  are  stored  in  the  array  WMACH(15) 
in  the  labeled  COMHOH  block  SOLMCH,  and  are  defiiierl  as  follows. 


WMACH(l) 

WMACH(2) 

WMACH(3) 

WMACH(4) 

WMACH(5) 

WMACH{6) 

WMACH(7) 

WMACH(8) 

WMACH(IO) 

WMACH(ll) 


is  NBASE,  the  base  of  floating-point  arithmetic. 

is  NDIGIT,  the  nundrer  of  NBASE  digits  of  precision. 

is  EPS,  the  floating-])oint  precision. 

is  RTEPS,  the  .s<iuare  root  of  EPSMCH. 

is  RMIN,  the  smallest  positive  floating-point  number. 

is  RTMIN,  the  square  root  of  RMIN. 

is  RMAX,  the  largest  positive  floating-point  number. 

is  RTMAX,  the  square  root  of  RMAX. 

is  NIN,  the  file  number  for  the  input  stream. 

is  NOUT,  the  file  number  for  the  output  stream. 


Within  routine  MCHPAR,  the  machine  constants  are  set  one  of  two  ways,  depending  upon  the 
value  of  the  logical  variable  HDWIRE,  which  is  set  in-hne. 

If  HDWIRE  is  .FALSE,  (the  value  set  for  the  distributed  copy  of  MCHPAR),  the  machine  constants 
are  computed  aiit.initiiticnlly  for  the  machine  being  used.  If  HDWIRE  is  .TRUE.,  machine  constants 
appropriate  for  the  IBM  3G0  Series  are  assigned  directly  to  the  elements  of  WMACH. 

Before  selecting  the  method  of  fissigning  the  machine  constants,  you  should  note  the  following. 
The  coTtqmtation  of  the  machine  constants  will  always  generate  a  single  arithmetic  underflow,  and 
hence  some  appropriate  remedial  action  may  need  to  be  taken  if  your  machine  traps  underflow. 

If  you  wish  to  implement  the  in-line  assignment  of  machine  constants  for  a  machine  other  than 
one  from  the  IBM  300/370  Series,  MCHPAR  mtist  be  modified  as  follows. 

1.  Change  the  in-line  assignment  of  HDWIRE  from  .FALSE,  to  .TRUE.. 

2.  S('t  the  values  of  WMACH  appropriate  for  the  machine  and  precision  being  used.  The  vfvlucs  of 
NBASE,  NDIGIT,  EPSMCH,  RMIN  and  RMAX  for  several  machines  are  given  in  the  following  table, 
for  both  single  and  double  precision;  RTEPS,  RTMIN  and  RTMAX  may  be  computed  using  Fortran 
statements.  The  values  NIN  and  NOUT  depend  on  the  machine  installation. 

For  each  precision,  we  give  two  values  for  EPSMCH,  RMIN  and  RMAX.  The  first  value  is  a  For¬ 
tran  decimal  approximation  of  the  exact  quantity;  use  of  this  value  in  MCHPAR  should  cause 
no  difficulty  except  in  extreme  circumstances.  The  second  value  is  the  exact  mathematical 
representation. 
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Table  of  machiiie>dcpendciit  parameters 


IBM  360/370 

Single 

CDC  6000/7000 

Single 

DEC  10/20 

Single 

— 

Univac  1100 

Single 

DEC  Vax 

Single 

NBASE 

16 

2 

2 

2 

2 

NDIGIT 

6 

48 

27 

27 

24 

EPS 

9.54E-7 

16'® 

7.11E-15 

2-47 

7.46E-9 

2-27 

1.50E-8 

2  26 

1.20E-7 

2-23 

RMIN 

l.OE-78 

16-65 

l.OE-293 

2-975 

l.OE-38 

2-129 

l.OE-38 

2-129 

l.OE-38 

2-128 

RMAX 

l.OE+75 

16®®(1-16'®) 

l.OE+322 

2l070(i_2-48) 

l.OE+38 

2l27(i.2-2T) 

l.OE+38 

2127(1-2-27) 

l.OE+38 
2127(1^2  2®) 

NBASE 


NDIGIT 


EPS 


IBM  360/370 

CDC  6000/7000 

DEC  10/20 

Double 

Double 

Double 

16 

2 

2 

14 

96 

62 

2.22D-16 

2.53D-29 

2.17D-19 

16->® 

2-95 

2-62 

l.OD-78 

l.OD-293 

l.OD-38 

16-85 

2-975 

2-129 

l.OD+75 

l.OD+322 

l.OD+38 

16®®(1-16*®) 

2l070(i_2-96) 

Univac  1100 
Double 


2 


61 


8.68D-19 

2-60 


l.OD-308 

2-1025 


l.OD+307 


DEC  Vax 
Double 


56 


2.78D-17 

2-55 


l.OD-38 

2-128 


l.OD+38 
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8.  EXAMPLE  PROBLEMS 

This  section  describes  a  linear  least-squares  problem  and  a  quadratic  program;  the  sample  main 
program  LSMAIN  that  calls  LSSOL  and  the  output  arc  given  in  the  Appendix. 

The  first  problem  is  a  constrained  least-squares  problem  of  type  LSI  with  nine  variables  and 
three  general  linear  constraints.  The  least-squares  matrix  and  vector  of  observations  arc  given  by 


/ 

1 

1 

1 

1 

1 

1 

1 

1 

1\ 

/1\ 

1 

2 

1 

1 

1 

1 

2 

0 

0 

1 

1 

1 

3 

1 

1 

1 

-1 

-1 

-3 

1 

1 

1 

1 

4 

1 

1 

1 

1 

1 

1 

1 

1 

1 

3 

1 

1 

1 

1 

1 

and 

b  = 

1 

1 

1 

2 

1 

1 

0 

0 

0 

-1 

1 

1 

1 

1 

1 

0 

1 

1 

1 

1 

1 

1 

1 

1 

0 

1 

1 

1 

1 

1 

1 

1 

1 

0 

1 

1 

1 

2 

2 

3 

1 

\ 

1 

0 

1 

1 

1 

1 

0 

2 

2  ' 

u/ 

The  least-squares  matrix  has  rank  6.  Let  t  in  LCLS  be  partitioned  into  two  sections:  the  first  n 
components  (denoted  by  tg),  corresponding  to  the  bound  constraints;  and  the  last  components 
(denoted  by  corresponding  to  the  linear  constraints.  The  vector  u  is  partitioned  in  a  similar 
fashion.  Using  this  notation,  the  upper  and  lower  bounds  on  the  variables  ate  given  by 


/b  —  (  2,  2,  oo, 

Ug  =  {  2,  2,  2, 


and  the  general  constraints  are  given  by 


The  starting  point  Xg  is 


,  -2,  -2,  -2,  -2,  -2)^ 

,  2,  2,  2,  2,  2f. 


1 

1 

1 

1 

^  oo  \ 

2 

1 

1 

1 

1 

1  and  «!,  = 

A 

1 

1 

1 

1 

l) 

[-2) 

xo  =  (  .1,  .5,  .3333,  .25,  .2,  .1667,  .1428,  .125,  .1111 


and  F{xo)  =  9.4746  (to  five  fig\ires).  The  optimal  solution  (to  five  figures)  is 

X  =  (  2.0000,  1.5719,  -1.4454,  -.037003,  .546685,  .17512,  -1.6567,  -.39477,  .31002  )\ 


and  F{x)  =  1.390587.  All  three  general  linear  constraints  are  satisfied  exactly  at  i*  The  Lagrange 
multiplier  associated  with  the  third  general  constraint  is  of  the  order  of  the  machine  precision,  and 
therefore  the  point  x  is  a  weak  minimum,  i.c.,  the  optimal  objective  function  is  unique,  but  is 
achieved  for  infinitely  many  values  of  x. 
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es 


The  second  problem  is  a  quadratic  programming  problem  of  type  QP2  with  a  semi-definite 
Hessian  matrix  and  linear  term  given  by 


(Note  that  by  setting  M  =  5,  we  need  not  assign  the  last  four  rows  and  columns  of  A  to  zero.) 

The  upper  and  lower  bounds  on  the  variables  are  given  by 

/„  =  (-2,  -2,  -2,  -2,  -2,  -2,  -2,  -2,  -2)^ 

lifl  =  (  2,  2,  2,  2,  2,  2,  2,  2,  2)'*'. 

and  the  general  constraints  are  given  by 

(-2\  /  1  1  1  1  1  1  1  1  4\  /1.5\ 

-2  ,  C  =  I  I  2  3  4  -2  1  1  1  1  I  and  =  1.5  I  . 

-2/  \  1  -1  1  -1  1  1  1  1  1  /  \  4  / 

The  starting  point  Xq  is  the  zero  vector,  at  which  F{xe)  —  0.  The  optimal  solution  (to  five  figures) 
is 

X*  =  (  2.0,  -.23333,  -.26667,  -.3,  -.1,  2.0,  2.0,  -1.7777,  -.45555  )^, 

and  /’(x*)  =  —8.067778.  The  first  two  linear  constraints  are  satisfied  exactly  at  the  solution,  as  are 
the  upper  bounds  on  variables  Xi,  xe  and  17.  Note  that,  although  the  Hessian  matrix  is  positive 
semi-definite,  the  point  x*  is  unique. 
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APPENDIX.  SAMPLE  PROGRAM  AND  OUTPUT 


Z  »  FILE  LSMAIN  FORTRAM 

3  * 

4  »  Sample  proaraa  for  Version  1.0  January  1964. 

5  +  -f-f +  ««  +  ++'» ««««««««« 

6 

7  IMPLICIT  DOUBLE  PRECISI0N(A-H,0-Z) 

8 

9  »  Set  the  declared  array  dimensions. 

10  *  NROMC  =  the  declared  row  dimension  of  C. 

11  •  KROWA  =  the  declared  row  dimension  of  A. 

12  »  MAXN  =  maximun  no.  of  variables  allowed  for. 

13  •  MAXM  =  maximum  no.  of  observations  allowed  for. 

14  *  MAXBND  =  maximum  no.  of  variables  +  linear  constraints. 

15  *  LIWORK  =  the  lerstth  of  the  integer  work  array. 

16  »  LMORK  =  tlie  leiigth  of  the  double  precision  work  array. 

17 

18  PARAMETER  (NROWC  =  3.  KROWA  =  10. 

19  $  MAXN  =  9,  MAXM  =  10, 

20  $  LIWORK  =  60.  LWORK  =  900, 

21  $  MAXBND  =  MAXN  +  NROMC  ) 

22 

23  INTEGER  KX(MAXN),  ISTATEI MAXBND 1 

24  INTEGER  IHORKI LIWORK ) 

25  DOUBLE  PRECISION  C( NROWC, MAXN),  BIMAXM) 

26  DOUBLE  PRECISION  BL(  MAXBND),  BUI  MAXBND),  CLAMOAtMA,XBM)) 

27  DOUBLE  PRECISION  CVEC(MAXN) 

23  DOUBLE  PRECISION  AIKROWA.MAXN) ,  X(MAXN) 

29  DOUBLE  PRECISION  KORKI LWORK) 

79 

3)  DOUBLE  PRECISION  BIGSND 

32  Ci:ARACTER»10  CBGSl® 

33 

34  INIRINSIC  FLOAT 

35 

36  PARAMETER  (  POINT) =0 . 1D*0,  POINT3=0.3D+0,  ONEPT5=1 .5D+0  ) 

37  PARAMETER  (  ZERO  =0.0D+0,  ONE  s1.00*0,  TWO  =2.00+0  ) 

33  PARAMETER  (  THREE  =3.00+0.  FOUR  =4.00+0,  FIVE  =5.00+0  ) 

39  PARAMETER  i  SIX  =6.00+0  ) 

40 

41  BIGEND  =  1.0D+15 

42  CCGBND  =  ' 1 .00+15* 

45 

44  •  ==r=r==i=i===r:====i=rs===r=i==z======r=== 

45  *  Example  1.  A  linear  least-squares  problem. 

46  »  ==:zz==r====z=====zz=zz=z=z===========r=== 

47  «  Set  tlie  actual  problem  dimensions. 

48  +  M  z  the  rRr;ber  of  observations  (rows  of  A)  (may  be  0). 

49  •  H  =  the  number  of  variables. 

50  *  hCLIN  =  the  number  of  general  linear  constraints  (may  be  0). 

51 

52  M  =10 

53  N  =9 

54  NCLIN  =  3 

55  NGM)  =  N  +  NCLIN 


in  DO  140  J  =  1,  N 

112  DO  130  1  =  1>  NCLIN 

113  C(I,J)  =  ONE 

114  130  CONTINUE 

115  140  CONTINUE 

116 

117  C(1,9)  =  FOUR 

113 

119  C(2,21  =  TWO 

120  C(2,3)  =  THREE 

121  C(2,4)  =  FOUR 

122  C(2,5)  =  -  TVS) 

123 

124  C(3,2)  =  -  ONE 

125  C(5,4)  =  -  ONE 

126 

127  DO  150  J  =  I,  N 

126  0L(J1  =  -  TWO 

129  BUtJi  =  TWO 

130  150  CONTINUE 

131  BL(  3)  =  -  BIGBND 

132 

133  «  Set  the  ranses  for  the  general  constraints. 

134 

135  BL(N+n  =  TWO 

136  EU(N»n  =  BIGBND 

137  EUN*2I  =  -  BIGBND 

138  BU(N+2)  =  -  TWO 

139  BL(N+3)  =  -  FOUR 

140  BU(N+3)  =  -  TWO 

141 

142  DO  170  J  =  1,  N 

143  X<J)  s  ONE  /  FLOAT(J) 

144  170  CONTINUE 

145 

146 

147  •  - 

148  *  Read  the  options  file. 

149  *  Add  a  single  option  using  a  call  to  LSOPTH. 

150  *  - 

151 

152  CALL  LSFILEI  lORTNSi  INFORM  1 

1D3  IF  (INFOT.M  .NE.  0)  THEN 

154  WRITE  (NOUT.  3000)  INFORM 

155  STO? 

156  END  IF 

157 

150  CALL  LSCPTNI  'Infinite  Bound  size  =V/CBG0ND  ) 

159 

163  *  - 

161  »  Solve  the  problew. 

162  •  - 

163 

164  CALL  L5SOL  (  M,  N, 

165  S  NCLIN,  NRCWC,  NROWA, 


so 
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166 

$ 

C>  BL>  BU,  CVEC, 

167 

$ 

ISTATE,  KX,  X,  A.  B, 

163 

t 

infcru,  iter,  obj,  claida. 

169 

$ 

IMORK,  LIU'ORK,  ts-ORK,  LUORK  I 

170 

171  • 

Test  for  an 

error  condition. 

1  72 

173 

IF  (INF(»M 

.GT.  1)  GO  TO  999 

Exsriple  Z.  A  QP  Mith  Hessian  bordered  by  zeros. 


Set  the  nsM  problem  dincnslons. 

M  =  the  nuniber  of  rous  (and  colunnsi  of  A  (may  be  0). 

N  =  the  number  of  variables. 

IXtlN  =  tite  nuo'.ber  of  scneral  linear  constraints  (may  be  0) 

CVEC  =  the  linear  part  of  the  objective  function. 


163 

184 

M  =5 

185 

N  =9 

156 

»CLIN  =  3 

167 

KBKO  =  N 

♦  NCLIN 

163 

189 

DO  220  J  = 

1,  M 

190 

DO  210  I 

=  1,  J- 

191 

A(I,J 

)  =  ONE 

192 

210 

CONTINUE 

193 

220 

COimSUE 

194 

195 

00  230  I  = 

1,  M 

196 

A(I,I)  = 

TWO 

197 

230 

CONTINUE 

198 

199 

DO  260  J  = 

1,  H 

200 

BL(J)  = 

-  TWO 

201 

BU(J»  s 

TWO 

202 

260 

CONTIh'UE 

203 

204 

BL(N*1  )  =  - 

TWO 

205 

B'J(N*1  1  = 

0NEPT5 

206 

BL(N»2)  =  - 

TWO 

207 

BU(N*2I  = 

0NEPT5 

208 

BL(N*3)  =  - 

TWO 

209 

BU(N*3)  = 

FOUR 

210 

211 

DO  270  J  = 

1,  N 

212 

CVEC(  J1 

=  -  ONE 

213 

270 

CONTINUE 

214 

CVEC( 1 )  =  - 

FOUR 

215 

CVEC(8)  =  - 

FOINTI 

216 

CVEC(9)  =  - 

P0INT3 

217 

218 

DO  280  J  = 

1,  N 

219 

X(J)  =  ZERO 

220 

280 

CONTINUE 

APPENDIX.  SAMPLE  PROGRAM  AND  OUTPUT 


226 

227  CALL  LSOPTNC  ‘Defaults  '  I 

228  CALL  LSOPTH(  ‘Problem  type  QP2‘  i 

229  CALL  LSOPTN(  ‘Rank  tolerance  =  I.OE-tO*  1 

230  CALL  LSOPTIK  ‘Feasibility  tolerance  a  t.0E-10‘  ) 

231 

232  *  - 

233  »  Solve  tSe  QP  problem. 


234  » 

235 

235 

CALL  LSSOL  ( 

M,  N, 

237 

$ 

NCLIN,  NROMC,  NROUA. 

233 

$ 

C.  BL.  BU>  CVEC. 

239 

$ 

ISTATE,  KX,  X,  A.  B, 

240 

$ 

INFORM,  ITER,  OBJ,  CLAMDA, 

241 

$ 

IimSK,  LIMORK,  MORK,  LNORK 

242 

243  * 

244 

Test  for  an 

error  condition. 

245 

IF  (INFORM  . 

GT.  1 )  GO  TO  999 

246 

STOP 

247 

2^3  »  Error  condition. 

249 

250  999  WPITE  (NOLtTi  30101  INFORM 

251  STOP 

252 

253  3000  FORHATi/  ‘  LSFILE  terminated  mith  INFORM  a*.  13) 

254  3010  FORMATi/  ‘  LSSOL  terminated  Mi th  INFORM  =*.  13) 

255 

256  *  End  of  the  example  program  for  LSSOL. 

257 

253  END 


ss 
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OPTIONS  flU 

BE6IN  Opilem  for  LSSOL  1.0  Sai^lo  prabloa. 

Iforoflons  LlolL  ■  £5 

Probloa  Typo  =  Laoat  oquaroo 

End 

Calls  fo  LSOPTN 

Infinlfo  Bound  sizo  *1.00*15 


SOL/LSSOL  -  Version  1.0  Fob  1086 

SS3=ZS==3SSS5SSSS33«SXS&SSSa8S3SSK&Sa 


Paranclers 


Problee  type . 

Linear  constraints . 

LSI 

3 

Feasibility  tolerance.. 

1.49E-0B 

COLO  start . 

Voriablcs . 

9 

Infinite  bound  size.... 

I.00E*I5 

Crash  tolerance . 

...  I.OOE-02 

objective  aatrix  roMS.. 

10 

Infinite  step  size..... 

1.D0E*I5 

Rank  tolerance . 

. . .  1 .49E-06 

EPS  leacblne  proclsioni 

£.22E-I6 

Feasibility  phase  Itns. 

60 

Print  level . 

...  to 

Korkspaeo  provided  la 

To  solve  problee  ee  need 

INI  60  li 

INI  9). 

Optlaallty  phase  Itns. 

Ml  900). 

Ml  261). 

25 

Rank  of  Iho  objcctlvo  function  data  oatrlM  ■  6 


I  In 

Jdel 

Jadd 

Step  NInf 

Slnf/Objectl ve 

find 

Lin 

Nz 

Nzl 

Nora  Gf 

More  Gzl 

Cond  T 

Cond  Rzl 

0 

0 

0 

0.0E*00 

2 

9.474603E*00 

0 

0 

9 

0 

6.86E*00 

0.00E*00 

l.0E*00 

0.0E*OO 

1 

1Z 

10L 

1.2E*00 

2 

5.987698E*00 

0 

1 

8 

0 

6.86E*00 

0.00E«00 

1.0E*00 

0.0E*0D 

2 

IZ 

nu 

4.IC-01 

1 

4.990079E*00 

0 

2 

7 

0 

3.00E*00 

O.00E*0O 

1 .1E*00 

0.0E*D0 

3 

1Z 

12U 

3.7E*00 

0 

4.95904IE*01 

0 

3 

6 

6 

5.60E*0I 

4.I3E*0) 

2.3E*00 

2.2E*0I 

4 

0 

tu 

3.0E-0I 

0 

2.429930E*0t 

1 

3 

5 

5 

3.89E*0I 

2.85E*0I 

2.4E*00 

4.8E*00 

5 

0 

0 

l.0E*00 

0 

1.390587E-0I 

1 

3 

5 

5 

6.55E-01 

1.59E-I5 

2.4E*00 

4.8E*0D 

Exit  froo  LS  problao  after  5  Iterations.  INFORN  *  1 


Variable 

State 

Value 

Lower  bound 

Upper  bound 

Lagr  aultipller 

Residual 

VAPBL 

f 

UL 

2.000009 

-2.000009 

2.000000 

-0.1191932 

0.0000E*00 

VARBL 

2 

FR 

1 .571959 

-2.000000 

2.000000 

o.oooooooE*oo 

0.4280 

VARBL 

3 

FR 

-1.445403 

None 

2.000000 

0.'"->q000E*00 

3.445 
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S3 


VARBL 

4 

FR 

-0.3700Z7SE-0I 

-Z. 000000 

Z. 000000 

O.OOOOOOOE^OO 

1.963 

VAPDL 

5 

FR 

0.5466054 

-Z. 000000 

Z. 000000 

O.OOOOOOOE^OO 

1.453 

VAPDL 

6 

FR 

O.I75tZ36 

-Z. 000000 

Z. 000000 

O.OOOOOOOE^OO 

1 .425 

VARBL 

7 

FR 

-1 .656704 

-Z. 000000 

Z. 000000 

O.OOOOOO0E«OO 

0.3433 

VA'fDL 

a 

FR 

-0.394774Z 

-Z. 000000 

Z. 000000 

o.oooooooE«oo 

1 .605 

VAROL 

9 

FR 

0.3I00Z90 

-Z. 000000 

Z. 000000 

o.oooooooE«oo 

1 .690 

Linear 

constr 

State 

Value 

Lomer  bound 

upper  bound 

Lagr  multiplier 

Residual 

L»)COH 

1 

LL 

Z. 000000 

Z. 000000 

None 

0. 39731 07E-0I 

-0.3553E-14 

UlCOtl 

Z 

UL 

-Z. 000000 

None 

-Z. 000000 

-0.1191932 

-0.42I9E-14 

L^^cor^ 

3 

UL 

-Z. 000000 

-4.000000 

-Z. 000000 

O.Z006660E-15 

-0.444IE-I5 

Exit  LSSOL  -  Weak  LS  solution. 

Final  LS  objoctlvo  value  »  0. 1390547 


C.ills  to  LSOPTN 


Defaults 

Prohlcn  type  QP2 

Pank  tolerance  =  I.OE-IO 

Feasibility  tolerance  ~  I.OE-IO 


SOL/LSSOL  -  Version  1.0  Feb  1904 


Parameters 


Problem  type . QPZ 

Linear  constraints .  3 

V..ri,.blcs .  9 

Cbjoctivc  eatrix  rotas..  5 

EPS  (inachine  precisioni  Z.2ZE-I6 


Feasibility  tolerance.. 
Infinite  botaad  size.... 
Infinite  step  size..... 


1.00E-10 

I.00E«I0 

I.OOE^IO 


COLD  start . 

Crash  tolerance........  1.00E-OZ 

Rank  tolerance . .  I.OOE-IO 


Feasibility  phase  Itns. 
Opt  leal I ty  phase  Itns. 


Print  level. 


Uorkspace  provided  Is  lUI 
To  solve  problem  we  need  1U( 


60)>  Ml  9001. 

9).  HI  Z70). 


Rank  of  the  objective  function  data  matrix  = 


Itn  Jdel  Jadd  Step  NInf  SInf/ObJectI ve  Bnd  Lin  Hz 


0 

0 

O.OE»00 

0 

O.OOOOOOE«00 

0 

tu 

7.5E-0I 

0 

-4.3750a0E«00 

0 

0 

1 .0E*O0 

0 

-4.400000E*00 

5? 

10U 

3.0E-0I 

0 

-4.700000E*00 

0 

0 

1 .OE»00 

0 

-5.I00000E*00 

7Z 

12U 

S.4E-01 

0 

-6.0S57I4E*00 

0 

6U 

1 .IE-02 

0 

-6.II3326E«00 

Norm  6f 
9.70E*00 
l.53E*00 
I .45E+00 
I .45E«00 
Z.47E*00 
Z.47E*00 

z.z-*oo 


Nora  Gzl 
4.47E+00 
5.00E-0I 
3.67E-17 
e.9<iE-0l 
I .ZOE-17 
I .73E+00 
I .64E«00 


Cond  T  Cond  Rzl 
Z.4E*00  t.3E«00 


Z.4E*00 
Z.4E+00 
l.0E«00 
I .0E«00 
Z.OE^OO 
Z.OE*00 


l.3E«00 
1 .3E*00 
1 .0E»00 
I .0E«00 
1 .3E«00 
I .7E*00 
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n 


7 

0 

IIU 

I.IE 

-01  0 

-6.ai5049E«00 

a 

3  4 

a 

a.o3E«oo  i.ieE*oo 

a.iE*oo 

1 . 

8 

0 

0 

1  .OE 

♦00  0 

-6.S30008E*00 

a 

3  4 

a 

i.ioE«oo  a.aaE-16 

a.iE*oo 

1 . 

9 

3Z 

0 

I.OE^OO  0 

-6.S67373E«00 

a 

3  4 

3 

1.07E«00  a.a3E-16 

a.iE*oo 

a. 

10 

9Z 

7U 

I.7E+00  0 

-8.0556iaE>00 

3 

3  3 

3 

3.83E-01  a.eCE-Ol 

a.iE*oo 

3. 

11 

0 

0 

l.OE^OO  0 

-0.067718E«00 

3 

3  3 

3 

4.3SE-01  1.05E-16 

a.iE«oo 

3. 

la 

IZU 

0 

I.OE^OO  0 

-6.067778E«00 

3 

a  4 

4 

4.3tE-01  t.OSE-16 

t.ZE^OO 

5. 

Exit  fro*  1 

OP  problaa  attar  IZ  Itaratlons.  INFORM  » 

0 

Variabla 

Stata 

Valua 

LoMar  bound 

Upper  bound 

Laar  aulti pilar  1 

Residual 

VARBL 

1 

UL 

a. 000000 

-a. 000000 

a. 000000 

-0.8000000 

O.OOOOE^OO 

VARDL 

a 

FR 

-0.a333333 

-a. 000000 

a. 000000 

O.OOOOOOOE^OO 

1  .767 

VARBL 

3 

FR 

-0.7666667 

-a. 000000 

a. 000000 

o.oooooooE«oo 

1.733 

VARBL 

4 

FR 

-0.3000000 

-a. 000000 

a. 000000 

o.oooooooE«oo 

1.700 

VARBL 

5 

FR 

-O.IOOOOOOE^OO  -a. 000000 

a. 000000 

O.OOOOOOOE^OO 

1  .900 

VARBL 

6 

UL 

a. 000000 

-a. 000000 

a. 000000 

-0.9000000 

0.0000E«OO 

VARBL 

7 

UL 

a. 000000 

-a. 000000 

a. 000000 

-0.9000000 

O.OOOOE^OO 

VARBL 

8 

FR 

-1.777778 

-a. 000000 

a. 000000 

o.oooooooE«oo 

o.aaaa 

VARBL 

9 

FR 

-0.4555556 

-a. 000000 

a. 000000 

O.OOOOOOOE^OO 

1.544 

Linear 

constr 

Stata 

Valua 

LoMar  bound 

Upper  bound 

Lagr  aultlplter 

Resichjal 

LUC  ON 

1 

UL 

1 .500000 

-a. 000000 

1.500000 

-0.6666667E-01 

-0.3553E-14 

LUCOH 

a 

UL 

1 .500000 

-a. 000000 

1.500000 

-0.33333S3E-01 

o.aaaoE-15 

LNCON 

3 

FR 

3.933333 

-a. 000000 

4.000000 

O.OOOOOOOE^OO 

0.6667E-01 

5E*00 

5C*00 

7E*00 

7E*00 

7E*00 

eE«oo 


Exit  LSSOL  -  Optiaal  QP  aolution. 

Final  QP  objactiv*  valua  x  -8.067778 


INDEX 


A  (ol)j('rtivc  data  luatrixjf  1. 

<*stiiiiatcd  rank  of,  4,  16,  18. 

iilontirally  zero,  16  (also  see  Linear  program). 

A.  9  (ilcfinition). 

Algoritliiu  of  LSSOL,  description,  2-6. 
n  (step  length),  2,  4 
printed  value,  17. 
u^f  (step  to  nciirest  constraint),  5. 

Amdahl  470,  20. 

ANSI  (1977)  Fortran,  1,  20. 

Artificial  constraint,  4  5  (definition),  17. 
Artificial  iimltiplier,  5. 

ASCII,  20. 

b  (vector  of  observations),  1. 

B.  9  10  (definition). 

Begin  (in  options  file),  12-13. 

BIGBND,  IS  (also  see  Infinite  Bound  Size). 

BL,  7  8  (definition). 

DLAS,  21. 

Level  2,  21. 

Bnd,  2,  17. 

BU.  8  (definition). 

Burroughs  6600  and  7600,  20. 

C  (general  constraint  matrix),  1. 

in  exiunples,  24  -25. 

Cm,  2,  3. 

Cm,  3. 

C.  7  (definition). 

CDC  6000  and  7000,  20. 

Checklist  of  optional  parameters,  16. 

Cholesky  factor,  3,  4,  9,  15. 

printout  cf  diagonals,  15. 

CLAMOA,  10  (definition). 

Cold  Start,  8,  14  (definition). 

Coliinui  interchanges,  4  (also  see  Rank). 
Coirunent  (in  optional  parameter  specification), 
12. 

Common  blocks,  list  of,  21. 

Cond  Rzl,  5,  18. 

Cond  T,  5,  18. 

Condition  estimator 
for  III,  5,  18. 
for  T,  5,  18. 

Condition  of  working  set,  control  of,  5-6. 
Constrained  stationary  point,  3. 

Con.str.iint  status  indicator  (see  ISTATE). 
Constr.-unt  violations,  weighted  sum  of,  17. 
Convexity,  2. 

Crash  Tolerance,  14  (definition). 

Cray-1,  20. 

CVEC,  8  (definition). 

Cyber,  20. 

Cycling,  10,  19. 

Data  General  MV/8000,  20. 

Data  matrix  (see  A  and  A). 

DEC  Systems  10  and  20,  20. 

DEC  VAX.  20. 

Default  values  of  optional  parameters,  checklist 
of,  16. 


Diagonals 

of  /?,  printout,  15. 
of  T,  |>riiitont,  15. 

Distribution  tape,  format  of.  20. 

DOUBLE,  7. 

Double  precision 

table  of  machine  constants,  23. 
version  of  code,  20. 

E  (priiite<l  constraint  designation),  17. 

End  (in  options  file),  12  13. 

EPS,  22  (also  see  t). 
t  (machine  precision),  14,  22.  . 

EQ  (printed  constraint  status),  18,  22. 

Equality  constraint,  1,  8,  17,  18. 

Error  correction  procedures,  19. 

Estimated  rank 

of  A,  4,  16,  18  (also  see  Rank  Tolerance), 
of  Ri,  16  (idso  sec  Rank  Tolerance). 

Example  1  (a  least-squares  problem),  24. 
Example  2  (a  quadratic  program),  25. 

Example  problems,  24-25. 

Explicit  lincikT  term  in  objective  function,  1,  10. 
External  file,  use  for  option  specification,  12-13. 

/  (transformed  residual),  3,  4. 

F  (objective  function),  1. 

Facom,  20. 

Feasibility  phase,  2,  4,  8,  17,  19. 

Feasibility  Phase  Iteration  Limit,  10,  14  (defi¬ 
nition). 

Feasibility  Tolerance,  2,  5,  8,  15  (definition), 
19. 

adjustment  to  avoid  overflow,  - 19. 

Feasible  point,  15  (definition). 

Feasible-point  problem,  1  (also  see  FP). 

Final  solution,  printout,  15.  ' 

Fixed  variables,  2  (also  see  EQ). 

Formal  parameters  of  LSSOL,  7-10. 

Formal  specification  of  LSSOL,  7. 

Format  of  distribution  tape,  20. 

Fortran 

ANSI  (1977),  20. 

subroutines,  naming  convention,  21. 

FP  (problem  type),  1,  7,  8,  16. 

FR  (printed  constraint  status),  18. 
fVee  variables,  2. 

Flyitsu,  20. 

Gabor,  Zsa  Zsa,  19. 

General  constraints,  1,  18. 

Global  minimum,  1. 

H  (Hessian  matrix),  1. 

BDWIRE,  22. 

Hessian  matrix,  1,  4. 

semi-definite  example,  25. 
upper-triangular  factor,  3,  4,  9. 

Hitachi,  20. 

Honeywell,  20. 

IBM 

360/370  and  3033/3081,  10,  20,  22. 

VS  Fortran,  20. 
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U:L  2900,  20. 

liiipU'inciitiition  infdrmation,  20  23. 
liifc'asibk'  i>r(>tilciii,  3  4,  10,  17,  19. 
Infoiusibilitip.s,  weighted  sum,  17. 
bilinite  lower  or  upper  bound,  1,  8. 

Infinite  Bound  Size,  8,  10,  15  (definition). 
Infinite  Step  Size,  15  (definition). 

INFORM,  10  (definition). 

Initial  working  set,  5,  6,  8  (also  see  Cold  Start 
and  Warm  Start). 

Input  parameter,  invalid,  10. 

Installation  procedure,  20. 

Invalid  input  par.amcter,  10. 
lOPTNS  (options  file  number),  12-13. 

ISTATE,  8  9  (definition),  14,  18. 

printout,  15. 

ITER,  10  (definition). 

Iteration  Limit,  15  (definition),  10. 

Iters  (see  Iteration  Limit). 

Itn  (printed  value),  17. 

Itns  (sec  Iteration  Limit). 

IW,  10  (definition). 

Jadd  (printed  value),  4,  17. 

Jdel  (i)rinted  value),  4,  17. 

Keyword  in  option  specification,  12. 

KX,  2,  9  (definition),  16. 

f  (lower  bound  vector),  1,  8  (also  see  BL). 

Ib,  24. 
t,,.  24. 

L  (printed  constraint  designation),  17  (also  see 
BL). 

Lagr  multiplier  (printe<i  value),  18. 

Lagrange  multiplier,  3,  10,  15,  18,  19,  24. 
optimal,  3  4,  10,  18. 
zero  or  small,  19. 

LCLS  (problem  statement),  1. 

Least  Squares  (sec  LSI). 

Le,a.st-s<|u.arrs  matrix,  1,  9  (also  see  A  and  A). 
Le.ast-squarcs  problem,  1. 

example.  24. 

LENIW,  11  (definition). 

LENW,  11  (definition). 

Level  2  DLAS,  21. 

Lin  (printed  value),  2,  17. 

Linear  constr,  18. 

Linear  least-sqtiiires  problem,  1. 

Linear  objective  function,  16. 

Linear  progr.am,  1,  16. 

Linear  Program  (sec  LP). 

Linear  term  in  objective  function,  1,  10. 

Lines  of  code  in  LSSOL,  1,  20. 

LL  (printed  constraint  status),  18. 

LNCON,  18. 

Local  minimum,  1  (also  sec  Weak  minimum). 
Lover  Bound,  17.  18  (also  see  BL). 

LP  (problem  type),  1,  7,  8,  16. 

LSI  (problem  type),  1,  7,  8,  16,  24. 

LS2  (problem  type),  1,  7,  8,  16. 

LS3  (problem  type),  1,  7,  8,  16. 

LS4  (problem  type),  1,  7,  8,  16. 

LSFILE,  12  13. 


LSOPTN,  13. 

list,  sample,  16. 

LSq  (.see  LSI). 

LSSOL 

idgorithm  of,  2  6. 
lines  of  code  in,  1,  20. 
parameters  of,  7-11. 
specification  of,  7. 

m,  1. 

mi,  (number  of  general  constraints),  1,  2,  5,  24. 
mw  (number  of  general  constraints  in  working 
set),  2. 

M,  7  (definition). 

Machine  constants 

computation  of,  21. 
tables  of,  22. 

Machine  dependencies  in  code,  21-23. 

Machine  precision  (see  c). 

Matrix  factorizations,  routines  for  updating,  21. 
MCHPAR,  22  (also  see  Machine  constants). 

Method  of  LSSOL,  de.scription,  2  -6. 

Minimal  sum  of  infeasibilities,  4,  17,  19. 
Minimum  abbreviation  (of  optional  parameter), 
14. 

n  (number  of  variables),  1. 
npR  (number  of  free  variables),  2. 
r»KX  (number  of  fixed  variables),  2. 
rin ,  3,  5,  17. 

N,  7  (definition). 

Naming  convention  for  Fortran  subroutines,  21. 
NDASE,  22. 

NCLIN,  7  (definition)  (also  see  mt). 

NDIGIT,  22. 

Negative  steps,  6  (also  see  a). 

NIN,  22. 

Nlnf  (number  of  infeasibilities),  17. 

No  feasible  point,  4,  10,  17,  19. 

Nolist  option,  13. 

Non-existent  lower  or  upper  bound,  8. 

None  (in  printout),  18. 

Nonlinearly  constrained  optimization,  6. 

NOUT,  22. 

Norm  Gf .  3.  17  (also  sec  Projected  gradient). 

Norm  Gzl,  5,  17  (also  see  Projected  gradient). 
NPSOL,  6. 

NROWA,  7  (definition). 

NROWC,  7  (definition). 

Null  space,  3. 

dimension  of  (sec  nz). 

Number  of  Lnfe.asibilitics,  17. 

Nz,  3,  5,  17. 

Nzl,  5,  17. 

OBJ.  10  (definition). 

Objective,  17. 

Objective  function  (F),  1. 
data  matrix  (see  A  and  A), 
linear,  16. 

Objective  matrix  (see  A  aud  A). 

Observation  vector  (5),  1. 

Optimal  Lagrange  multiplier,  3-4  (definition), 
10,  18. 


Optimal  solution,  10. 

Optimality  pluase,  2  (also  see  MetlioJ  of  LSSOL). 
Optimality  Phase  Iteration  Limit,  1(1,  14  (defi¬ 
nition). 

Optimality  test,  10. 

O]>tion-handling  routines,  21. 

Optional  parameters,  12-16. 
checklist  and  default  values,  10. 
cumulative  changes,  13. 
description,  14-16. 

Oi;tions  file,  12  -13. 

Ordering  of  variables,  2  (also  see  KX). 

Orthogonal  transformation,  2. 

Overflow,  19. 

p  (search  direction),  2,  3. 

Pkr,  2,  3. 

Pariuneter  vector  (see  x). 

Parjuneters  of  LSSOL,  7-11. 

Phase  1  (see  Feasibility  phase). 

Pluise  2  (see  Optimality  phase). 

Phrase  (to  modify  optional  parameter),  12. 
Precision,  machine  (see  e). 

Primal  method,  2. 

Prime  Systems,  20. 

Print  Level,  10.  11,  IS  (definition). 

Printed  output,  de.scription,  17-18. 

Printout,  control  of,  15. 

Problem  type  (see  Problem  Type). 

Problem  Type,  1,  15  -16  (definition). 

Projected  gradient,  3,  10. 
norm,  17. 

Q,  3 
Qfr,  2. 

QP  (problem  type),  7,  8. 

QPl  (]>roblem  type),  1,  7,  8,  16. 

QP2  (problem  type),  1,  7,  8,  16,  25. 

QP3  (problem  type),  1,  7,  8,  16. 
qP4  (problem  type),  1,  7,  8,  16. 

QR  factorization,  4. 

Quadratic  program,  1,  1C. 
example,  25. 

Quadratic  Program  (problem  type)  (see  qP2). 
Qu.alifying  phrase  (in  optional  parameter),  12. 

R,  3,  4,  9,  15  (also  see  iZi). 
ordering  of  columns  (see  KX). 
printout  of  diagonals,  15. 

Ri,  4,  16. 

condition  estimate  of,  18  (also  see  Cond  Rzl). 
Rz,  3,  4,  10,  17. 

singular,  3,  4,  10,  17. 

Rank,  4,  16. 

determination,  16. 

Rank  Tolerance,!,  10  (definition),  18. 

REAL,  7. 

References,  26. 

Re-ordering  of  v.ariables,  2  (also  see  KX). 

Reset  optional  par.aineters,  13-14. 

Residual,  10,  18. 

Residual  vector,  3,  4,  10. 

Reverse-triangular  matrix,  2  (also  see  T). 

RHAX,  22. 


RMIN,  22. 

RTEPS,  22. 

RTMAX,  22.- 
RTNIN,  22. 

Search  direction  (p),  2. 

Second-derivative  matrix,  1  (also  see  Hcssi.an  ma¬ 
trix). 

Semi-definite  Hessian  matrix,  example,  25. 
Sequential  quadratic  ]>rogramming  method,  6. 
Simplex  method,  2,  5. 

Simplex  steps,  2,  5. 

Sinf  (weighted  sum  of  infeasibilities),  17. 

Single  precision 

table  of  machine  constants,  23. 
version  of  code,  20. 

Singuhir  Rz,  3,  4,  10,  17. 

Small  Lagrange  multiplier,  10. 

Source  files,  list,  20. 

Specification  of  LSSOL,  7. 

Standard  simplex  method,  2. 

State,  18  (also  see  ISTATE). 

Stationary  point,  3. 

Status  intlicator  fur  constraints  (see  ISTATE). 

Step  (juinted  value),  17  (also  sec  Step  length). 
Step  length  (a),  2,  4,  17. 
choice  of,  4,  5. 

Strong  local  minimum,  1,  10. 

Sum  of  infeasibilities,  3-4. 
minimum,  17. 
weighted,  17. 

Synonyms  (for  optional  parameters),  12. 

T,  2,  5. 

condition  estimate,  18  (.also  see  Cond  T), 
printout  of  di.agunals,  15. 

T  (printed  constraint  designation),  17  (also  see 
Temporary  bound). 

Tape 

characteristics,  20. 
format  20. 

Temporary  bound,  5,  17. 

TQ  factorization,  2,  5,  15. 

TVansformed  residual  vector  (/),  3,  4,  10. 
IVapezoidal  matrix,  1,  9  (.also  sec  TViangular  fac¬ 
tor). 

Triangular  factor,  3,  4,  9,  15. 

of  Hessian  as  data  matrix,  0. 

Two-phase  )>rinial  method,  2. 

u  (upper  bound  vector),  1,  17  (adso  see  BO). 
u«,  24. 

UI..24. 

0  (printed  constraint  designation),  17. 

UL  (printed  constraint  status),  18. 

Unbounded 

objective  function,  10,  15. 
solution,  1,  10. 
step,  15. 

Underflow,  19. 

Unique  solution,  1,  10. 

Uuivac  1100,  20. 

Unknowns,  vector  of  (see  z  and  X). 

Updating  matrix  factorizations,  routines  for,  21. 


Upper  bound,  18  (also  sec  u  and  BU). 
Uppor-trapezoidiU  matrix,  1,  9  (also  sec  'IViang- 
ular  fiM-tor). 

Upper-triangnlar  factor  (see  Triangular  factor). 

Valid  option  strings,  examples  of,  12. 

Value,  18. 

VARBL,  18. 

Variable, 18. 

Variance-covariance  matrix,  9. 

Vector 

of  observations  (6),  1  (also  sec  B). 
of  unknowns  (z),  1  (also  see  X). 

Vertex,  2,  3,  17, 

Violations,  constraint  (see  Infeasibilities). 

W,  11  (definitinn). 

Warm  Start,  6,  8,  14  (definition),  19. 

Weak  minimum,  1,  10. 
example  of,  24. 

Weak  LP  solution  (see  Weak  minimum). 

Weak  LS  solution  (see  Weak  minimum). 

Weak  QP  solution  (sec  We.ak  minimum). 
Weighted  sum  of  infeiisibililics,  4,  17  (also  see 
bifcasible  problem). 

WMiCH,  22  (also  see  Machine  constants). 

Working  precision,  7  (also  sec  e). 

Working  set 
changes  in,  4. 

condition  estimate,  18  (also  see  Cond  T). 
definition,  2. 

Workspace  parameters  of  LSSOL,  10-11. 

z  (vector  of  unknowns),  1. 
printout,  18. 

X,  0  (definition). 

y,  2. 

Z  (basis  for  null  space),  2  (also  sec  Null  space). 

dimension  of  (see  nz). 

Zi,  4-5,  17. 

^iSrR’  5  (also  sec  Projected  gradient). 

Zj,  4-5. 

Z  (printed  constraint  designation),  17. 

Zero  Lagrange  multiplier,  19  (also  see  Lagrange 
multiplier). 

--  (printed  constraint  status),  18  (also  see  Infea¬ 
sible  problem). 

♦  +  (printed  constraint  status),  18  (also  sec  Infea¬ 
sible  problem). 
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